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ABSTRACT 


This  report  describes  the  development  of  an  advanced  multitarget  tracking 
algorithm  for  surveillance  radars.  The  approach  is  based  on  the  development 
of  mathematical  models  relating  targets,  the  environment,  and  radar  parameters 
to  the  measurements  supplied  to  the  tracking  algorithm.  The  fundamentals  of 
the  surveillance  and  tracking  algorithm  are  based  on  adaptive  hypothesis  test¬ 
ing  techniques,  and  specifically,  on  the  multiple  model  adaptive  estimation 
approach.  The  algorithm  was  developed,  computer-programmed,  and  exercised  on 
several  hypothetical  but  realistic  surveillance/ tracking  scenarios. 
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SECTION  1 
OVERVIEW 

This  report  is  concerned  with  the  development  and  evaluation  of  an  ad¬ 
vanced  tracking  algorithm  for  surveillance  radars.  Similar  algorithms  have 
been  studied  for  many  years.  Recent  advances  in  computer  technology  make  a 
re-examination  of  current  tracking  algorithms  timely,  since  it  will  be  possi¬ 
ble  to  implement  more  sophisticated  algorithms  in  the  near  future.  The  most 
promising  algorithms  are  based  on  statistical  estimation  theoretic  concepts, 
specifically  on  adaptive  hypothesis  testing  techniques.  Most  of  the  individ¬ 
ual  Issues  that  arise  in  the  design  of  tracking  algorithms  (false  and  missing 
returns,  target  maneuvers,  data  association,  etc.)  have  been  previously 
addressed  within  this  framework,  but  research  is  needed  to  handle  all  of  these 
issues  simultaneously. 

The  algorithm  developed  in  the  present  research  effort  provides  the  capa¬ 
bility  to  incorporate  most  of  the  subjective  conditions  associated  with  a 
realistic  surveillance/tracking  scenario.  At  present,  the  algorithm  can 
handle  a  multitarget  tracking  situation,  it  can  identify  targets’  "births”  or 
"deaths”  (i.e.,  appearance  or  disappearance)  and  distinguish  these  from  a 
detection  or  a  missed  detection.  In  the  future,  the  algorithm  might  be  aug¬ 
mented  with  additional  features.  It  could  detect  and  estimate  the  occurrence 
and  magnitude  of  target  maneuvers,  and  it  could  have  the  capability  of  identi¬ 
fying  target  formations,  missile  launches,  etc. 
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The  development  of  the  tracking  algorithm  Is  firmly  grounded  in  statis¬ 
tical  estimation  theory;  it  is  an  example  of  the  multiple  model  adaptive  esti¬ 
mation  approach.  In  its  optimal  form,  the  algorithm  cannot  be  implemented  due 
to  data  storage  and  processing  requirements  which  grow  without  bound.  But 
starting  from  the  complete  optimal  algorithm,  we  developed  rules  for  limiting 
these  requirements  without  significantly  compromising  tracking  performance  and 
implemented  the  resulting  algorithm.  Using  a  selection  of  simple  but  realis¬ 
tic  scenarios  the  power  of  the  tracking  algorithm  and  its  resiliency  in  the 
face  of  stressing  environments  has  been  demonstrated.  The  implementation  is 
quite  general  and  flexible,  providing  an  excellent  example  of  how  to  organize 
a  multitarget  algorithm  for  stressing  radar  surveillance  environments. 

This  report  begins  with  the  development  of  a  mathematical  model  of  the 
target  state  equations  and  of  the  radar  measurements  (including  such  environ¬ 
mental  factors  as  clutter  and  multipath).  It  should  be  emphasized  that  this 
model  is  not  intended  to  be  a  so-called  "truth  model"  capable  of  faithfully 
replicating  the  precise  details  of  radar  returns,  but  rather  a  "design  model" 
used  to  ensure  that  all  of  the  environmental  issues  are  reflected  in  the 
tracking  algorithm.  Subsequently,  we  develop  the  so-called  algorithmic  or 
simulation  model.  This  model  describes  the  target/radar  relationship  in  a 
manner  suitable  for  synthetic  radar  measurement  generation  which  are  needed  to 
exercise  the  surveillance/tracking  algorithm.  These  models  are  developed  and 
described  in  detail  in  Section  2. 

Once  the  radar  environment  models  are  established,  they  are  incorporated 
in  the  surveillance  and  tracking  algorithm.  The  optimal  version  of  the  track¬ 
ing  algorithm  is  devel sped  in  Section  3;  it  constitutes  a  member  of  the  so- 
called  adaptive  estimation  algorithms  family.  Section  3  begins  with  a 
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description  of  the  most  basic  technique  for  adaptive  estimation,  the  multiple 
model  estimation  algorithm.  The  optimal  algorithms  are  next  presented  for  the 
three  cases  of  environment  phenomena.  These  cases  have  the  following 
formulation : 

1.  one  target  with  measurement  errors,  missed  detections,  and 
false  detections,  but  no  maneuvers; 

2.  one  target  with  measurement  errors,  missed  detections, 
false  detections,  appearance  and  disappearance,  but  no 
maneuvers;  and 

3.  an  unknown  number  of  targets  with  measurement  errors, 
missed  detections,  false  detections,  appearance  and  dis¬ 
appearance,  but  no  maneuvers. 

Each  algorithm  is  more  complex  than  the  preceding,  but  the  progression  through 
the  cases  aids  understanding  of  the  algorithms.  Section  3  finishes  with  a 
useful  reformulation  of  the  algorithm  for  the  most  complex  case;  it  is  the 
formulation  which  is  used  in  implementing  the  algorithm. 

At  this  point  it  will  have  become  apparent  that  implementing  the  optimal 
full-blown  tracking  algorithm  is  infeasible.  This  is  evident  from  the  expo¬ 
nentially  growing  number  of  data  association  hypotheses  with  each  additional 
scan.  Clearly,  practicality  dictates  that  the  number  of  hypotheses  actually 
generated  and  tested  by  a  tracking  algorithm  be  bounded  by  some  value  depen¬ 
dent  on  the  host  computer.  Consequently,  we  develop  in  Section  4  methods  that 
relax  the  optimal  algorithm.  The  relaxation  rules  include  several  approaches 
such  as  ignoring  unlikely  hypotheses  (e.g.,  a  "birth"  immediately  followed  by 
a  "death"),  gating  of  data  associations  (this  rule  "gates-out”  returns  from 
consecutive  scans  which  are  unreasonably  separated,  thus  implying  impossible 
target  measurement  errors),  and  pruning  (i.e.,  only  a  fixed  number  of  the  most 
likely  hypotheses  is  allowed). 
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In  Section  5  the  computer  implementation  of  the  algorithm  is  discussed. 
In  Section  6  the  computer  program  is  exercised  on  four  hypothetical  (but 
realistic)  scenarios  involving  two  and  three  crossing  targets.  It  is  then 
demonstrated  that  the  (theoretically)  suboptimal  algorithm  (i.e.,  relaxed  by 
the  rules  developed  in  Section  4)  performs  as  accurately  as  the  optimal  ver¬ 
sion.  Finally,  conclusions  drawn  from  this  research  effort  are  summarized  in 
Section  7. 
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SECTION  2 

MODELING  OF  TARGETS  AND  OF  THE  RADAR  ENVIRONMENT 

Research  into  tracking  theory  for  airborne  surveillance  radars  requires 
two  sorts  of  radar  environment  models.  One  sort  is  mathematical  and  describes 
the  relationship  of  the  radar  returns  to  target  variables  (position,  velocity, 
etc.)  as  a  function  of  the  radar  environment  (thermal  noise,  clutter,  etc.)  in 
a  manner  suitable  for  tracking  algorithm  design.  The  other  sort  of  model  is 
algorithmic  and  describes  the  same  relationship  in  a  manner  suitable  for  syn¬ 
thetic  radar  measurement  generation.  Synthetic  radar  measurements  are  needed 
to  exercise  any  tracking  algorithm  and  to  determine  its  performance  and 
robustness.  In  this  section,  we  first  present  our  algorithm  design  model. 

Then  we  briefly  present  the  very  simple  radar  measurement  generation  model 
used  to  produce  the  algorithm  performance  reported  in  Section  6.  F.'nally,  we 
sketch  a  design  for  a  more  sophisticated  radar  model  which,  if  implemented, 
would  allow  more  thorough  testing  of  airborne  surveillance  radar  tracking 
algorithms . 

2.1  ALGORITHM  DESIGN  MODEL 

Several  target-relative  phenomena  should  be  modeled:  the  number  of  tar¬ 
gets  (aircraft  and  missiles)  within  radar  coverage  can  vary  from  scan  to  scan 
as  aircraft  land,  take  off,  crash,  are  shot  down,  or  fly  into  or  out  of  cover¬ 
age;  in  effect  targets  are  "mortal,''  being  "born”  and  “dying."  Targets  can 
move  on  essentially  constant  headings  with  essentially  constant  velocities  or 
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they  can  maneuver  (within  reasonable  limits).  An  aircraft  can  launch  a  mis¬ 
sile,  i.e.,  one  target  splits  into  two,  or  be  hit  by  one,  i.e.,  two  targets 
merge  into  one.  In  addition,  several  radar-related  phenomena  should  be 
modeled:  each  target  can  produce  only  one  detection  per  scan,  but  may  produce 
none  because  of  thermal  noise,  clutter,  scintillation,  etc.  False  detections 
may  be  caused  by  thermal  noise,  clutter,  jamming,  etc.  Errors  occur  in  mea¬ 
suring  target  locations  due  to  thermal  noise,  clutter  jamming,  glint,  etc. 
Targets  flying  in  formation  may  be  unresolvable . 

Each  phenomena  included  requires  a  model.  Some  models  can  be  quite  sim¬ 
ple,  as  in  the  measurement  errors  which  can  be  treated  as  white  Gaussian  noise 
of  known  covariance.  Others  can  be  relatively  complex,  as  in  target  maneu¬ 
vers.  And  many  possible  models  are  plausible  for  some  phenomena,  such  as  tar¬ 
get  "mortality,"  i.e.,  appearance  and  disappearance.  The  intuitive  model  of 
target  "mortality"  is  simple;  but  a  formal  model  is  needed  specifying  the 
probabilities  of  various  numbers  of  targets  being  "born"  per  scan,  the  proba¬ 
bility  density  function  of  target  states  at  "birth,”  the  probabilities  of  tar¬ 
gets  "dying"  on  a  given  scan,  etc. 

Our  model  development  effort  has  proceeded  by  considering  just  a  few 
phenomena  (without  trivializing  the  problem) ,  then  adding  additional  phenomena 
one  step  at  a  time.  The  cases  we  have  formulated  are: 

1.  one  target  (known)  with  measurement  errors,  missed  detec¬ 
tions,  and  false  detections,  but  no  maneuvers  or  other 
complex  phenomena; 

2.  one  target  (known)  with  measurement  errors,  missed  detec¬ 
tions,  false  detections,  "mortality,”  but  no  maneuvers, 
etc.;  and 

3.  an  unknown  number  of  targets  with  measurement  errors, 
missed  detections,  false  detections,  "mortality,"  but  no 
maneuvers,  etc. 
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Complete  models  for  these  three  cases  have  been  developed  and  are  described  in 
the  succeeding  three  subsections.  Presentation  of  the  models  in  this  manner 
makes  the  third  (and  most  complex)  far  more  comprehensible  than  it  would  be  if 
it  were  presented  in  isolation. 

2.1.1  Case  1 

The  first  requirement  is  for  a  model  of  the  target's  dynamic  state.  We 

use 

x(t)  *  Ax(t-l)  +  w(t-l)  for  all  t  (2-1) 

where  w(t)  is  Gaussian  distributed  with  mean  0  and  variance  W.  We  can  further 
assume  that  the  components  of  the  dynamic  state  vector  are  target  position 
North  and  East  of  some  arbitrary  point  and  the  North  and  East  target  veloci¬ 
ties.  Then  we  would  have 

x(t)  «  lP£(t) 

and 

A  - 


W  - 


pnU)  PeU)  PnU)]t 


(2-2a) 


WaTs 


(2-2b) 
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where  I2  Is  the  2x2  Identity  matrix,  Ts  is  the  scan  time,  and  Wa  is  the  covar¬ 
iance  of  the  East  and  North  target  accelerations  produced  by  wind,  etc.  Note 
that  it  is  possible  to  define  the  algorithm  without  reference  to  the  dimension 
or  content  of  the  dynamic  state  vector  or  of  the  covariance  of  the  state  driv¬ 
ing  noise. 

The  dynamic  state  of  the  target  is  not  its  complete  state.  At  any  given 
moment,  the  one  target  may  or  may  not  be  detected;  if  detected,  it  may  be  any 
one  of  the  detections.  We  define  the  target's  detection  state  as  d(t),  which 
consists  of  two  components: 


d(t)  -  [b(t),i(t)]T  (2-3) 

where 

b(t)  e  {D,D}  (2-4) 

and 

i(t)  e  {],... ,L( t) }  (2-3) 

where  b(t)  equals  D  if  the  target  is  detected  at  t  and  D  otherwise,  where  i(t) 
equals  the  index  of  the  target  detection  if  the  target  is  detected  at  t,  and 
where  L(t)  is  the  number  of  detections  at  t. 

The  detection  state,  like  the  dynamic  state,  is  Markovian  because  we 
assume 

Prob[d(t) |d(t)  for  T<t]  “  Prob[d(t)J  for  all  t  .  (2-6) 

The  basic  property  of  conditional  expectations  allows  us  to  write 
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Probld(t)]  *  Prob( i(t) |b(t) ]Prob[b(t) )  for  all  t  ,  (2-7) 

and  we  model  the  two  terms  on  the  right-hand  side  of  Eq.  2-7  as  follows: 


Prob[b(t)]  -  1 

for 

b(t)-D 

and 

L(t)-0 

1-PD 

for 

b(t)»D 

and 

L(t )>0 

for 

b(t)-D 

and 

L(t)>0 

0  otherwise  , 


and 


Prob[ i(t) |b(t)]  -  1/L( t )  for  b(t)«D  and  l<i(t)<L(t) 
0  otherwise  . 


(2-9) 


The  complete  state  of  the  target  is  the  hybrid  of  the  continuous-valued  dynam¬ 
ic  state  x  and  the  discrete-valued  detection  state  d.  This  hybrid  is  also 
Markovian  because  we  assume  the  dynamic  and  detection  states  evolve 
independently. 

The  single  target  is  assumed  to  cause  at  most  one  detection  per  scan,  but 
any  number  of  false  detections  may  occur.  We  denote  the  number  of  false  de¬ 
tections  at  t  by  Lp^(t)  and  assume  it  to  be  Poisson  distributed  with  parameter 
Xp^.  That  is, 

-*FA  LFA(t) 

Prob[LpA(t)]  -e  Ufa  /LfaU)1!  for  LfaU)>0 


0 


otherwise 


(2-10) 
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L(t)  -  LFA(t)  if  b(t)  -  D 


(2-11) 


L(t)  -  LFA(t)+l  otherwise, 


so  L(t)  is  conditionally  Poisson  distributed.  Also  note  that  XpA  is  the 
expected  number  of  false  alarms  per  scan. 

The  above  paragraphs  describe  the  detection  process  for  the  1  target  in 
the  presence  of  false  alarms,  but  not  the  measurement  process.  We  define 
z(i,t)  to  be  the  measurement  of  the  i-th  detection  at  t  and  let 


z(i(t),t)  -  Cx(t)+v(t)  if  b(t)=D  (2-12a) 

z(k,t)  -  v'(k,t)  for  ke{l , . . . ,L(t ) }-{ i( t) }  if  b(t)“D  (2-12b) 

z(k,t)  ■  v'(k,t)  for  ke{l , . . . ,L( t)}  if  b(t)*D  (2-12c) 

where  v(t)  is  Gaussian  distributed  with  mean  0  and  covariance  V  while  v'(i,t) 
is  uniformly  distributed  over  the  surveillance  set  S.  This  surveillance  set 
is  a  closed  subset  of  the  measurement  space.  We  further  assume  that 


C  -  [i2  0] 


(2-13) 


and 


E  m  (z  |  Izl  <  maximum  surveillance  range}  .  (2-14) 

Again,  it  Is  not  necessary  to  be  so  restrictive  in  defining  the  algorithm;  we 
have  chosen  this  simple  model  for  our  initial  development  in  order  to  avoid 
the  obscuring  of  the  key  concepts  to  be  developed. 
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It  will  prove  convenient  to  let  |s|  symbolize  the  area  or  volume  of  the 
surveillance  set.  Then 


Prob[v'(i.t)]  -  l H | — 1  if  v’(i,t)  e  = 
0  otherwise 


(2-15) 


where  Prob[v ' (i  ,t )  ]  is  shorthand  for  the  probability  density  function  of 
v'(i,t).  It  will  also  prove  convenient  to  define  a  constant,  pp^i  such  that 

PFA  *  *FA/ |=|  •  (2-16) 


This  constant  is  the  expected  number  of  false  detections  per  scan  and  per  unit 
area  or  volume  and  will  be  referred  to  as  the  false  alarm  density;  it  is  pro¬ 
portional  to  the  radar  false  alarm  rate. 

2.1.2  Case  2 

Taking  into  account  target  "mortality,"  we  must  model  the  "birth"  and 
"death"  of  the  single  target.  We  do  so  by  expanding  the  definitions  of  d(t) 
and  b(t)  in  the  previous  case.  As  before 

d(t)  -  [b(t),i(t)jT  (2-17) 

but  in  this  case, 

b(t)  e  {B,D,D,X}  (2-18) 

where  b(t)  equals  B  if  the  target  has  never  been  detected  on  or  before  t 
(i.e.,  it  is  "unborn"),  equals  D  if  the  target  is  detected,  equals  D  if  the 
target  is  undetected  but  has  been  detected  and  will  be  again  (i.e.,  it  is 
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still  “alive"),  and  X  if  the  target  will  never  be  detected  again  (i.e.,  it  is 
"dead").  These  definitions  of  the  values  which  b(t)  may  take  hint  at  the 
restrictions  on  the  state  transitions  allowed;  the  single  target  is  "born" 
only  when  it  is  first  detected  and  cannot  be  detected  again  once  it  dies. 

Because  it  now  involves  more  than  just  detection,  we  rename  d(t)  the 
event  state  of  the  1  target;  it  includes  all  of  the  discrete  components  of  the 
complete  state  of  the  1  target.  The  event  state  is  Markovian  because  we 
assume 

Prob[d(t) | d ( t )  for  r<t-l]  =  Prob[  i(t) |b(t) ]Prob[b(t) |b(t-l) J  for  all  t>0 

(2-19) 

Prob[ i(t) | b(t) ]  -  1/L(t)  for  b(t)=D  and  l<i(t)<L(t) 


(2-20) 


0  otherwise  , 


Prob[ b( t) I b(t-l) ] 


1-PB 

for 

b(t)-B 

and 

b( t-1 )“B 

Pb 

for 

b( t)“D 

and 

b(t-l)*=B 

Pd 

for 

b(t)=D 

and 

b( t-1 )=D 

i-Pd-px 

for 

b(t)«D 

and 

b(t-l)-D 

Px 

for 

b(t)-X 

and 

b(t-l)-D 

Pd 

for 

b(t)-D 

and 

b(t-l)-D 

1-PD 

for 

b(t)“D 

and 

b(t-l)-D 

1 

for 

b(t)-X 

and 

b(t-l)-X 

otherwise 


(2-21a) 


and 


12 


ALPHATECH,  INC 


— 

Problb(O)]  -  1  for  b(0)-B 
0  otherwise 

The  state  transitions  of  b(t)  are  diagrammed  in  Fig.  2-1. 


( 2—2  lb) 


Figure  2-1.  State  Transition  Diagram  for  b(t)  with  Associated 
Transition  Probabilities. 

The  figure  clearly  illustrates  that  we  have  intimately  coupled  "mortal¬ 
ity,"  a  target-related  phenomenon,  and  detection,  a  radar-related  phenomenon. 
On  the  surface,  it  would  seem  preferable  to  have  less  coupling;  to  associate  a 
target's  "birth"  not  with  its  first  detection  but  rather  with  its  entry  into 
the  radar's  coverage  (or  some  other  such  event  external  to  the  radar).  To  do 
so  would  cause  additional  complexity  in  the  resulting  tracking  algorithm. 
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Detection  is  an  observable  event  to  the  radar  tracking  algorithm,  but  target 
entry  into  the  radar's  coverage  is  not  observable  in  the  absence  of  target 
detection;  the  radar  receives  no  information  about  target  positions  if  the 
target  is  not  detected.  If  one  were  to  model  without  the  coupling  introduced 
here,  the  resulting  algorithm  would  take  into  account  the  possibility  of  ex¬ 
isting  targets  which  have  not  yet  been  detected.  Clearly  such  hypotheses  are 
indistinguishable  from  ones  which  avoid  hypothesizing  undetected  targets. 

This  point  illustrates  the  care  which  must  be  taken  in  modeling  the  radar  en¬ 
vironment;  models  which  seem  natural  or  intuitively  pleasing  on  the  surface 
may  introduce  subtle  but  substantial  complexities  in  the  tracking  algorithm. 

The  model  of  the  single  target's  dynamic  state  must  also  be  modified  for 
this  case.  It  becomes  the  following: 


Ax(t-l)+w(t-l)  if 

b(t)*D  or  D 

and 

b(t-l)«D  or  D 

xQ  if 

b(t)*D 

and 

b(t-l)«B 

(2-22) 

where  Xq  is  Guassian  distributed  with  mean  0  and  covariance  ZQ.  Effectively, 
xQ  is  the  dynamic  state  with  which  the  target  is  "born.”  Once  "born,"  the 
single  targets  dynamic  state  evolves  in  the  same  manner  as  Case  1  until  the 
target  "dies."  The  dynamic  state  of  the  target  before  its  "birth"  or  after 
its  "death"  is  undefined. 

2.1.3  Case  3 

In  this  case,  the  number  of  targets  is  greater  than  one.  We  model  the 
total  number  of  targets  as  indefinitely  large,  but  assume  that  a  finite  number 
are  "born"  on  any  given  scan.  Thus,  the  number  "born"  at  or  before  t,  denoted 
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by  N(t),  is  always  finite;  it  also  is  monotonically  increasing  with  time.  As 
with  the  previous  case,  it  is  necessary  to  extend  the  definition  of  d(t).  In 
this  case,  we  let 

d(t)  -  tN(t),d'(l,t)T,d*(2,t)T,...,d*(N(t),t)TlT  (2-23) 

where 

d'U,t)  *  lbU,t),\U,t))T  ,  (2-24) 

bU,t)  e  {B.D.D.X}  ,  (2-25) 

and 

\(i,t)  c  ,L(t)J  .  (2-26) 

Effectively,  d'(£,t)  is  the  event  state  of  the  i-th  target  where  that  state  is 
as  defined  for  the  I  target  in  Case  2.  The  total  discrete  valued  state  of  the 
system  of  targets  d(t),  is  finite  dimensional  because  it  only  includes  the 
Individual  event  states  d*(i,t)  of  the  N(t)  targets  "born*’  on  or  before  scan 
t,  but  the  dimension  grows  with  time. 

This  event  state  is  Markovian.  In  particular,  we  assume 

N(t) 

Prob(d(t)  |d(t-l))  -  Prob(N(t)|N(t-l))  II  ProblbU.t)  |b(i,t-l)  ,N(t)  ,N(t-l)  1 

t-1 

x  p(i(l,t) . \(N(t),t)jb(I,t) . b(N(t),t))  (2-27) 

where 
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Problb(i,t)|bU,t-l),N(t),N(t-l)] 


1 

for 

bU,t)«B 

and 

N(t)<* 

1 

for 

bU,t)-D 

and 

N(t-l)<fc<N(t) 

Pd 

for 

bU,t)«D 

and 

b(£,t-l)*=D 

i-Pd-Px 

for 

bU,t)-D 

and 

b(£,t-l)=D 

Px 

for 

b(l,t)=X 

and 

b( 4 ,  t-1 )»D 

Pd 

for 

b(£,t)=D 

and 

b(i,t-l)=D 

1-PD 

for 

bU,t)-D 

and 

b(i,t-l)=D 

1 

for 

b(i,t)=X 

and 

b( i ,t-l )=X 

0 

otherwise 

(2-28) 

,  ~AB  [ N( t )-N( t-1 ) J 

Prob[N(t) | NC t— 1 ) ]  *  e  XB  /(N(t)-N(t-l) ] !  for  N(t)>N(t-l) 

and  t>0  , 

( 2-29a) 


and 


Prob[N(0) ]  -  1  for  N(0)=0 

0  otherwise  (2-29b) 

We  will  define  the  last  term  in  Eq.  2-27  in  a  moment,  but  it  is  useful  to 
pause  at  this  point  to  provide  an  interpretation  of  the  above  expressions 
which  actually  are  simpler  than  they  may  appear.  The  state  transitions  of 
b(k,t)  are  diagrammed  in  Fig.  2-2.  This  diagram  is  virtually  Identical  to 
Fig.  2-1;  the  difference  is  that  the  transition  from  B  to  D,  i.e.,  the  "birth" 
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of  Che  target  is  governed  by  the  N(t)  process.  When  a  target's  number  comes 
up,  it  is  "born."  The  number  of  targets  "born"  on  a  given  scan  is  Poisson 
distributed  with  mean  Xg  and  is  independent  from  6can  to  scan. 


1  IF  N(t)<£ 


Figure  2-2.  State  Transition  Diagram  for  b(£,t)  with  Associated 
Transition  Probabilities. 

We  now  define  the  probabilistic  description  of  the  data  association  vari¬ 
ables  1  ( £ ,t) •  The  key  here  is  to  note  that  specifying  the  b(£,t),  £-1 . 

N(t)  determines  which  i(£,t)  must  be  defined  (i.e.,  which  targets  are  to  have 
measurements  associated  with  them  at  time  t).  It  is  useful  to  divide  this  set 
of  targets  into  three  groups: 
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IdbU)  - 

{£|b(£,t)-D&b(£,t-l)«B} 

(2-30a) 

*DD<t)  “ 

{£|b(£,t)«D&b(£,t-l)-D} 

( 2-30b) 

*DD  “  (*! 

|b(£,t)-D&b(£,t-l)«D}  . 

(2-30c) 

We  use  the  lower-case  letter  i^g(t),  and  to  denote  the 

cardinalities  of  these  sets.  The  probability  distribution  for  the  i(£,t)  for 
£  in  the  union  of  these  three  sets  can  be  described  as  follows:  of  the  total 
set  of  L(t)  returns  at  time  t,  it  is  equally  likely  that  any  subset  of 
cardinality 

>fT)B(  c  )  +  toDU)  +  (2-31) 

corresponds  to  the  set  of  returns  from  targets;  out  of  this  subset,  it  is 
equally  likely  that  any  subset  of  cardinality 

tDD(c)  +  <toD(t)  (2-32) 

corresponds  to  the  set  of  returns  from  targets  which  have  been  detected  at 
previous  times,  and  any  permutation  of  the  elements  of  this  subset  corresponds 
to  an  equally  likely  association  of  measurements  in  this  subset  with  the  set 
of  previously  existing  targets  which  are  detected  at  time  t.  The  remaining 
♦db(c)  returns  to  be  associated  correspond  to  returns  from  new  targets.  Note 
that  here  there  are  a  number  of  absolutely  indistinguishable  possibilities, 
corresponding  to  the  ordering  of  the  new  targets  in  our  list.  We  remove  this 
ambiguity,  as  we  describe  in  more  detail  later,  by  ordering  these  targets  in  a 
unique  manner,  corresponding  to  the  magnitude  of  the  range  measurement 
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associated  with  each.  Given  the  above  discussion,  one  can  compute  the  proba¬ 
bility  for  any  allowable  set  of  values  of 

i(£,t),ie'fDI(t)U’fDD(t)U'i'DD(t)  .  (2-33) 

Here  "allowable"  means,  for  example  that  no  two  of  the  i(A,t)'s  can  take  on 
the  same  value,  as  individual  returns  correspond  to  individual  targets.  The 
resulting  probability  is: 

p( i  ( 1 ,t),... , i (N( t ) ,t) | b ( 1 , t ) , . . . ,b(N(t),t)) 

-  p(  iU,t),UH'Di(t)UfDD(t)U'}'DD(t)) 


[L(t)-»DD(t)-¥DD(t)-fDB(t)] !¥DB! 


L(t)! 


(2-34) 


This  probability  is,  as  expected,  uniform  over  the  set  of  allowable  values. 

Just  as  it  was  necessary  to  define  a  total  discrete  valued  state  for  the 
system  of  targets,  made  up  of  individual  event  states,  it  is  necessary  to  de¬ 
fine  a  total  continuous  valued  state,  made  up  of  the  individual  target's  dy¬ 
namic  states.  Thus  we  let 


x ( t )  =  [N(t),x'(l,t)T,x'(2,t)T,...,x,(N(t),t)T]T  for  t>0  (2-35) 


where 

x'U.t)  *=  Ax'  U,t-l)+w(£,t-l)  if  b(*,t)-D  or  D  and  b(£,t-l)=D  or  D  , 

x'U)  if  bU,t)=D  and  b(Jt,t-l)=B  (2-36) 

o 
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and  where  x'(£)  is  independently  Gaussian  distributed  for  all  t  with  mean  0 
o 

and  covariance  I'.  Each  target's  dynamic  state  is  modeled  in  this  case  as  the 
o 

single  target's  dynamic  state  was  modeled  in  Case  2. 

The  measurement  model  must  also  be  modified  to  account  for  multiple  tar¬ 
gets  and  to  be  consistent  with  the  data  associations  specified  by  the  i(£,t) 
and  b(£,t).  Specifically, 

z(i(£,t),t)  -  Cx'U.t))  +  v(k,t)  for  i£fDB(t)U¥DD(t)UHlDD(t)  (2-37a) 

z(j  ,t)  =  v'(j,t)  for  jet(t)  ( 2—37  b ) 

where  v(k,t)  is  independently  Gaussian  distributed  with  mean  0  and  covariance 
V,  where  v'(j,t)  is  uniformly  distributed  over  the  surveillance  set  E,  as  in 
Case  1,  and  where 

4>(t)  -  { 1 , . . .  »L(t) }  -  (iU,t)  Ue1'DB(t)UH'DD(t)UTDD(t)}  .  (2-38) 

As  discussed  earlier,  in  order  to  remove  ambiguity  in  ordering  new  targets  we 
place  the  following  constraint  on  the  i(£,t):  should  contain  more  than 

one  element,  we  restrict  the  indices  it  contains,  N(t-l)  through  N(t),  to  be 
ordered  such  that 

lz(iU,t),t)l  <  lz(  iU’  ,t),t)»  (2-39) 

for  all  £  and  £'  in  f Dg( t )  such  that  £<£’.  That  is  to  say,  if  more  than  one 
target  is  "born"  at  t,  we  assume  they  are  indexed  in  magnitude  order  of  the 
measurements  they  produce.  Because  each  measurement  is  a  random  value,  this 
ordering  is  unique. 
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2.1.4  Environment  Model  Accuracy 

Certain  aspects  of  these  models  are  somewhat  unrealistic.  For  example, 
it  is  assumed  that  the  airborne  surveillance  radar  provides  the  tracker  with  a 
360°  "snapshot"  of  detected  target  positions  once  per  scan.  In  typical  air¬ 
borne  surveillance  radars,  each  azimuth  is  revisited  once  per  scan,  but  in  a 
rotating  search  pattern.  More  significantly,  position  measurement  accuracy, 
detection  probability,  and  the  distribution  false  detections  are  all 
assumed  to  be  independent  of  range  and  azimuth.  In  typical  airborne  surveil¬ 
lance  radars,  measurement  accuracy  degrades  with  range,  because  of  the  signal- 
to-thermal  noise  power  ratio  decreases  with  range  and  the  width  of  the  radar 
beam  (in  meters,  not  radians)  increases  with  range.  Also,  clutter,  and  the 
signal  processing  techniques  used  to  suppress  its  impact ,  can  have  a  very  non- 
uniform  effect,  blocking  target  detection  or  producing  dense  clusters  of  false 
detections  in  local  regions  of  the  surveillance  volume. 

The  obvious  question  is  whether  a  tracking  algorithm  based  on  a  somewhat 
unrealistic  model  can  perform  well.  There  is  reason  to  believe  that  an  un¬ 
realistic  model  need  not  produce  poor  tracking  performance.  In  simpler, 
single  target  tracking  problems,  Kalman  filters  are  often  applied  even  when 
the  uncertainties  in  target  motion  or  measurement  are  known  to  be  non-Gaussian 
in  distribution  and  of  variance  different  from  that  assumed.  Yet  these  Kalman 
filters  produce  position  and  velocity  estimates  and  predictions  which  are 
quite  satisfactory;  they  perform  quite  robustly.  There  may  even  be  an  advan¬ 
tage  to  basing  a  tracking  algorithm  on  a  somewhat  unrealistic  model  if  it 
reduces  the  sensitivity  of  model  performance  to  mismatches  between  the  actual 
and  modeled  environments.  Consider  a  tracking  algorithm  based  on  an 
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exceedingly  accurate  model  of  sea  clutter  returns  over  water;  its  performance 
may  be  better,  but  how  will  its  performance  be  degraded  by  ground  clutter  when 
used  over  land,  or  over  an  archipelago? 

The  question  of  whether  satisfactory  performance  can  be  achieved  by  a 
tracking  algorithm  based  on  a  somewhat  unrealistic  model  of  the  radar  environ¬ 
ment  cannot  be  definitively  answered  using  analytic  techniques.  The  use  of 
analytic  techniques  requires  the  assumption  of  an  environment  model  (possibly 
different  from  that  upon  which  the  tracking  algorithm  is  based),  leading  to 
the  question  of  how  accurate  this  model  must  be.  A  definitive  answer  requires 
the  use  of  actual  radar  data  to  drive  the  algorithm  and  "truth"  data  which 
reliably  reports  the  actual  positions  and  velocities  of  the  tracked  aircraft. 
Of  course,  such  data  is  expensive  to  obtain.  Simulated  data  can  be  clearly  as 
valuable.  Simulation  also  requires  use  of  an  environment  model,  but  allows  a 
very  complex  model.  By  using  a  significantly  more  complex  model  to  produce 
simulated  radar  data  than  the  model  upon  which  the  tracking  algorithm  is 
based,  the  robustness  of  the  tracking  algorithm  in  the  face  of  environment 
modeling  inaccuracies  can  be  tested. 

2.2  SIMULATION  MODEL 

As  noted  immediately  above,  an  environment  model  is  needed  in  any  radar 
simulation  producing  data  to  drive  a  tracking  algorithm  and  that  model  should 
be  more  complex  than  (or  at  least  different  from)  the  model  upon  which  the 
algorithm  is  based  in  order  to  test  the  algorithm's  robustness.  Early  in  this 
project  we  produced  a  high-level  design  for  an  ideal  airborne  surveillance 
radar  simulation  for  testing  tracking  algorithms.  In  concert  with  our  track¬ 
ing  algorithm  software  development,  we  began  developing  software  for  a 
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simplified  simulation  with  a  goal  of  building  toward  the  ideal  simulation. 

The  high-level  design  of  the  ideal  simulation  is  described  in  this 
subsection. 

Figure  2-3  illustrates  our  conceptual  model  of  the  context  in  which  the 
tracking  algorithm  operates.  We  have  decomposed  the  conceptual  model  into 
several  linked  submodels  which  can  be  usefully  isolated  intellectually.  The 
"external  signal  sources"  are  target  aircraft,  jammers,  and  ground  clutter; 
the  "external  signal  sources"  submodel  describes  these  sources,  including  tar¬ 
get  aircraft  motion.  Motion  of  the  airborne  early  warning  (AEW)  aircraft  is 
described  in  another  submodel.  All  pretracking  functions  and  phenomena  of  the 
AEW  radar  system  are  described  in  the  "radar  transformation  process"  and  the 
"radar  measurement  process"  submodels.  The  division  of  the  radar  system  into 
two  submodels  is  somewhat  arbitrary,  but  convenient.  Functionally,  the  divi¬ 
sion  is  between  the  output  of  the  Doppler  filters  and  the  input  to  the  detec¬ 
tion  processor.  In  order  to  suppress  ground  clutter  returns  in  the  antenna’s 
main  lobe,  the  "radar  measurement  process"  must  have  available  an  estimate  of 
the  AEW  aircraft  ground  speed  and  heading.  We  assume  that  estimate  is  pro¬ 
vided  by  the  "AEW  aircraft  navigation  system"  which  we  describe  rather  simply. 
The  figure  names  the  information  provided  by  each  submodel  to  the  others;  this 
information,  and  the  submodels,  are  briefly  described  in  the  next  few  para¬ 
graphs.  Time  is  treated  as  discrete,  not  continuous,  in  all  of  the  submodels. 

2.2.1  "External  Signal  Sources"  Submodel 

We  model  individual  target  aircraft  motion  as  arbitrary  within  realistic 
constraints  on  accelerations  and  velocities.  This  approach  proved  successful 
in  recent  ALPHATECH  research  for  MICOM  on  aircraft  maneuver  detection.  The 
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Figure  2-3.  Conceptual  Simulation  Model. 
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number  of  target  aircraft  is  also  considered  arbitrary  but  constrained  to  be 
reasonable.  Targets  flying  in  formation  can  be  modeled;  since  they  would 
typically  be  unresolvable ,  they  are  considered  collocated.  Jammers  may  also 
be  collocated  with  aircraft  or  they  may  be  stationary.  Ground  clutter  is 
modeled  as  both  diffuse  (producing  noise-like  returns)  and  discrete  (producing 
target-like  returns).  At  each  discrete  time  step,  the  "external  signal 
sources"  submodel  supplies  maps  of  the  signal  sources  describing: 

1.  current  target  aircraft  positions, 

2.  current  jammer  positions  and  effective  radiated  power, 

3.  diffuse  clutter  densities,  and 

4.  positions  of  discrete  clutter  sources. 

The  map  coordinates  are  "absolute"  in  that  they  are  rectangular  and  are  refer¬ 
enced  to  a  fixed  location,  not  to  the  current  position  of  the  AEW  aircraft. 

We  initially  are  making  several  simplifying  assumptions  about  jamming  and 
clutter.  Jammers  are  assumed  to  be  stationary,  all  of  equal  power,  and  uni¬ 
formly  distributed.  We  similarly  assume  that  discrete  clutter  sources  are  all 
of  equal  size  and  uniformly  distributed.  Diffuse  clutter  is  assumed  to  be  of 
uniform  density  everywhere. 

2.2.2  "AEW  Aircraft"  Submodels 

Motion  of  the  AEW  aircraft  is  assumed  to  be  periodic,  along  a  racetrack 
course.  The  "AEW  aircraft  motion"  submodel  also  outputs  a  map  in  "absolute" 
coordinates  at  each  discrete  time  step,  giving  the  current  AEW  aircraft  posi¬ 
tion,  orientation,  and  velocity.  The  "AEW  aircraft  navigation  system"  sub¬ 
model  presently  is  assumed  to  pass  these  values  through  without  alteration  at 
each  time  step.  In  reality,  there  are  errors  introduced  into  the  tracking 
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process  by  use  of  estimates  of  A£W  aircraft  position,  orientation,  and  veloc¬ 
ity.  But  we  believe  that  such  errors  are  small  compared  to  others  and  so  have 
chosen  not  to  model  this  source  at  this  time. 

2.2.3  "Radar  Transformation  Process"  Submodel 

The  AEW  aircraft's  radar  is  modeled  as  doubly  transforming  the  maps  of 
the  external  signal  sources  in  the  course  of  producing  the  Doppler  filter  out¬ 
puts.  The  first  transformation  is  a  geometric  one:  using  the  current  AEW 
aircraft  position,  the  maps  are  translated  into  radar-centered  polar  coor¬ 
dinates.  The  second  transformation  is  into  sufficient  statistics  at  the  out¬ 
put  of  the  Doppler  filters.  Actually,  the  outputs  of  the  Doppler  filters  are 
random  voltages.  But  we  are  not  interested  in  modeling  these  voltages  them¬ 
selves;  we  are  interested  in  modeling  the  detection  and  location  of  target 
aircraft  using  these  voltages.  Thus,  it  is  efficient  to  compute  only  statis¬ 
tics  sufficient  to  describe  the  random  voltages  from  the  viewpoint  of  the 
detection  and  location  processes.  Both  the  overall  model  and  the  simulation 
based  upon  it  are  simplified  by  using  sufficient  statistics  at  this  stage.  In 
the  "radar  transformation  process"  submodel  we  are  also  making  some  simplify¬ 
ing  assumptions  initially.  For  example,  all  antenna  and  Doppler  filter  side- 
lobes  assume  constant  magnitude  while  the  main  lobes  are  "squared  off." 

The  sufficient  statistics  are  described  by  two  range/azimuth/Doppler  fre¬ 
quency  maps.  One  is  for  target  returns  (including  multipath  returns)  and  dis¬ 
crete  clutter  returns.  It  gives  the  expected  range,  azimuth,  and  Doppler  of 
the  peak  return  as  well  as  the  expected  magnitude  of  the  peak.  The  second  map 
is  for  noise-like  interference  (jamming,  diffuse  clutter  returns,  and  receiver 
noise  combined).  It  gives  the  expected  interference  power  as  a  function  of 
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range,  azimuth,  and  Doppler.  The  sufficient  statistics  are  calculated  for 
each  azimuth  as  if  that  azimuth  were  at  the  center  of  the  beam  and  the  beam 
were  stationary.  Thus,  the  rotation  of  the  radar  antenna  and  all  its  effects 
are  not  modeled.  This  omission  should  have  little  impact  on  tracking  algo¬ 
rithm  performance  and  greatly  simplifies  the  model  (and  the  simulation). 

2.2.4  "Radar  Measurement  Process”  Submodel 

The  "radar  measurement  process"  models  all  detection  and  post-detection 
radar  system  features  short  of  the  tracker.  The  detection  algorithm  is  quite 
simple:  response  peaks  in  range,  azimuth,  and  Doppler  above  a  fixed  threshold 

are  declared  detections  except  those  occurring  in  the  azimuth/Doppler  band 
predicted  for  main  lobe  ground  clutter  returns  or  those  occurring  along  azi¬ 
muths  where  excessive  jamming  creates  a  very  high  interference  level.  Figure 
2-4  illustrates  the  azimuth/Doppler  band.  For  each  azimuth  value,  that  band 
is  constant  width  in  Doppler  and  is  centered  on  the  currently  estimated  AEW 
aircraft  ground  speed  along  that  azimuth.  Detection  statistics  can  be  calcu¬ 
lated  from  the  two  maps  of  sufficient  statistics  using  standard  formulae  for 
the  detection  of  signals  in  noise  or  noise-like  interference.  False  alarm 
statistics  can  be  similarly  calculated  from  the  map  of  sufficient  statistics 
for  noise-like  interference. 

Interpolation  of  detection  range  and  azimuth  is  modeled  by  assuming  inde¬ 
pendent  Gaussian  measurement  errors  with  variances  inversely  proportional  to 
mean  target  response-to-mean  interference  ratio.  The  resultant  polar,  AEW- 
relative  measurements  are  translated  into  "absolute"  coordinates  for  input  to 
the  tracker  using  the  current  estimate  of  the  AEW  aircraft  position. 
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Figure  2-4.  Azimuth/Doppler  Band. 

It  is  necessary  also  to  translate  the  measurement  covariance  into  abso¬ 
lute  coordinates.  Figure  2-5  illustrates  the  translation.  Azimuthal  uncer¬ 
tainty  must  be  translated  into  cross  range  uncertainty  and  then  the  down  range 
and  cross  range  uncertainties  must  be  translated  into  absolute  coordinate 
uncertainties . 

2.2.5  "Tracker"  Submodel 

The  tracking  algorithm  requires  estimates  of  target  aircraft  detection 
probabilities  and  of  false  alarm  rates  in  order  to  correctly  account  for 
missed  and  false  detections.  In  some  AEW  radar  systems,  these  quantities  may 
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Figure  2-5.  Measurement  Covariance  Transformation. 
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be  estimated  from  signal  processor  outputs  by  sophisticated  means  while  in 
other  systems  nominal  values  are  assumed.  For  simplicity,  we  have  chosen  to 
assume  nominal  values.  The  values  are  modeled  as  provided  by  the  "radar  mea¬ 
surement  process"  as  they  would  be  in  more  sophisticated  cases. 
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SECTION  3 

ALGORITHM  DEVELOPMENT :  OPTIMAL  ALGORITHM 

Optimal  linear  estimation  algorithms  are  by  now  well  understood  from  both 
a  theoretical  and  applied  point  of  view.  In  contrast,  aside  from  the  adhoc 
extended  Kalman  filter,  no  general  nonlinear  estimation  algorithms  are  avail¬ 
able.  But  for  a  specific  class  of  problems,  specifically,  problems  that  would 
be  linear  except  for  the  presence  of  certain  unknown  discrete  or  continuous 
parameters,  the  family  of  so-called  adaptive  estimation  algorithms  are  avail¬ 
able.  The  problems  of  target  tracking  in  the  environment  described  in  the 
previous  section  fall  into  this  class.  Such  algorithms  are  designed  to  be  op¬ 
timal  by  some  probabilistic  criterion,  e.g.,  Bayesian,  maximum  likelihood,  or 
maximum  a  posteriori  probability. 

This  section  begins  with  a  description  of  the  most  basic  technique  for 
adaptive  estimation,  the  multiple  model  estimation  algorithm.  This  basic 
technique  assumes  a  finite  number  of  unknown  discrete  parameters;  its  general¬ 
ization  to  a  number  of  discrete  parameters  which  grow  without  bound  is  then 
described.  The  optimal  algorithms  are  next  presented  for  the  three  cases  of 
environment  phenomena  in  the  previous  section  (see  pp.  5  and  6).  Each  algo¬ 
rithm  is  more  complex  than  the  preceding,  but  the  progression  through  the 
cases  aids  understanding  of  the  algorithms.  The  section  finishes  with  a  use¬ 
ful  reformulation  of  the  algorithm  for  the  most  complex  case;  it  is  the  formu¬ 
lation  which  was  used  in  implementing  the  algorithm. 
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3.1  MULTIPLE  MODEL  ADAPTIVE  ESTIMATION  (MMAE)  ALGORITHM 

3.1.1  Basic  Algorithm 

The  basic  MMAE  algorithm  is  defined  for  a  linear  system  of  the  form 
State  Dynamics 


x  (  t+ 1 )  =  A(  t  )  34  ( t  )  +  w(  t ) 

Measurement  Equation 

£(t+l)  =  C(t+l)34(t+l)  +  v(  t+1 ) 

where  t=0,l,2,.. 

.  is  the  time  index. 

x(  t  leR-1 

is  the  state  vector  (nonwhite  stochastic 
sequence ) , 

w(c)eRJ 

is  the  white  plant  noise. 

v(t)eRK 

is  the  white  measurement  noise,  and 

z(t)eRK 

is  the  measurement  vector. 

Probabilistic  Information 

The  initial  state  x(0)  is  Gaussian  with 

E{x(0)}»x(0)  (3-3) 

cov[x(0) ;x(0) ]  «  *  l'Q  >  0  (3-4) 

The  plant  noise  w(t)  is  Gaussian  discrete  white-noise  with 

E{w(t)}  -  0  (3-5) 

cov[w(t) ,w( t) ]  -  W(t)6tx  (3-6) 
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W(t)  ■=  WT(t)  >  0  .  (3-7) 

The  measurement  noise  _v(t)  is  Gaussian  discrete  white-noise  with 

E{v(t)}=0  (3-8) 

cov[v(t);v(i)]  =  V.(t)6tT  (3-9) 

V(t)  -  VT(t)  >0  (3-10) 

(i.e.,  every  measurement  is  corrupted  by  white  noise) 

x(0)  ,  w(t)  ,  _v(t)  (3-11) 

are  independent  for  all  t,t. 

It  is  assumed  that  the  system  of  concern  has  been  modeled  in  the  above 
form,  but  that  some  parameters  of  the  model  are  undetermined.  If  we  denote 
the  unknown  parameters  by  a  vector  _d,  the  dependence  of  the  model  matrices  can 
be  made  explicit  by  the  notation  A(t;d_),  W(t;d_),  etc.  it  is  assumed  that  the 

parameter  or  model  set  is  finite  and  discrete,  jd e *— 2 . — ' n)*  This  means, 

in  essence,  that  N  hypotheses  are  possible  concerning  the  system  and  that  the 
goal  is  to  identify  which  is  correct  as  well  as  to  estimate  the  system  state. 

The  solution  the  MMAE  estimation  algorithm  provides  is  a  means  to  find 
the  optimal  (conditional  expectation)  estimate  x(t|t)  of  the  state  x(t)  and  a 
means  to  compute  the  probability  p^(t)  that  the  i-th  parameter  is  the  correct 
parameter  (alternatively,  that  the  i-th  hypothesis  is  correct),  given  the  ob¬ 
servations  z(0) , . . .  ,z( t ) .  Several  authors  have  considered  this  problem;  the 
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original  derivation  is  apparently  due  to  Magill  [1].*  At  any  rate,  straight¬ 
forward  Bayesian  manipulations  show  that  the  solution  is  given  by  the 
relations 


Bi(t)exp 


-  rT(t)S-1(t)ri(t) 
2  ~ i  ~i 


PiU)  = 


) 


I  Bj(t)exp 

j  =  l 


-  -  £j(t)Ti_1(t)£j(t)  j  pj(t-l) 


) 


Pi(t-l)  (3-12) 


N  - 

x(t|t)  =  Z  x^(t\t)pi(t)  .  (3-13) 

i=  1 

In  Eqs.  3-12  and  3-13,  the  subscript  i  denotes  a  quantity  computed  assuming 
that  model  i  is  the  correct  model,  i.e.,  that  d= .  Thus  Xi.(t|t)  is  the 
standard  Kalman  filter  estimate  for  the  model  of  Eqs.  3-1  through  3-11  with 
d*=Dj ,  £i(t)  is  the  corresponding  filter  residual 

£i.  ( t )  *  zU)  -  C(t;Di)x^(t  1 1  —  1 )  ,  (3-14) 

_Si(t)  is  the  residual  covariance  matrix,  and 

Bi(t)  =  |SjL(t)|l/2  .  (3-15) 

The  solution  is  illustrated  in  Fig.  3-1. 


♦References  are  indicated  by  numbers  in  square  brackets  and  appear  at  the  end 
of  this  report. 
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Figure  3-1.  Multiple  Model  Adaptive  F.stlmation  Algorithm. 
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The  MMAE  algorithm  has  been  the  object  of  considerable  study,  both  analy¬ 
tic  and  by  simulation,  e.g.,  [2]  —  [7  J .  It  has  been  known  for  some  time  that  if 
the  true  system  is  in  the  model  set  and  is,  say,  model  k,  then  under  reason¬ 
able  assumptions 


Pk<t)  -  1  (3-16) 

Pi(t)  +0  ,  i  *  k  (3-17) 

as  t-*-c°.  More  interesting  is  the  realistic  situation  in  which  the  true  model 
is  not  in  the  model  set,  but  may  be  some  high  order,  even  nonlinear  system. 

In  this  case,  it  has  been  shown  that 

Pk(t)  -  1  (3-18) 

Pi(t)  +  0  ,  i  *  k  (3-19) 

where  k  minimizes  a  certain  information  distance  between  the  set  of  admissible 
models  and  the  true  system,  [8]-[10].  This  information  distance  is  readily 
computed  (in  the  case  in  which  the  true  system  is  linear)  by  standard  covari¬ 
ance  calculations  associated  with  mismatched  Kalman  filter  performance  evalua¬ 
tion.  Moreover,  it  can  be  shown  that  the  convergence  in  Eqs.  3-18  and  3-19 
(in  an  appropriate  probabilistic  sense  made  precise  in  [8])  is  exponential 
with  coefficient  equal  to  the  information  distance. 

3.1.2  Generalizations  of  MMAE  Algorithm 

Two  generalizations  of  the  above  described  MMAE  algorithm  are  necessary 
before  it  can  be  applied  to  multitarget  tracking.  The  first  of  these  restric¬ 
tions  is  that  the  dimension  J  of  the  state  vector  is  assumed  not  dependent  on 
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d  in  the  preceding  subsection.  But  there  are  many  situations  in  which  the 
dimension  of  the  state  vector,  i.e.,  the  order  of  the  system,  is  one  of  the 
unknown  parameters;  it  will  be  demonstrated  below  that  tracking  multiple  tar¬ 
gets  is  one  such  situation.  It  is  then  impossible  to  combine  the  conditional 
state  estimates  into  an  unconditional  estimate  as  per  Eq.  3-13.  Under  such 
circumstances,  a  maximum  a  posterior  probability  criterion  is  adopted,  i.e., 
for  each  t 

x ( 1 1 1 )  -  Xi(t|t)  (3-20) 

where  i  is  such  that 

Pl(t)  >  pj(t)  (3-21) 

for  all  j.  As  above,  if  the  true  system  is  in  the  model  set  and  is,  say, 
model  k,  then  under  similar  reasonable  assumptions, 

Pk(t)  -  1  ( 3-22a ) 

and 

Pj(t)  -*-0  for  i  *  k  (3-22b) 

as  t+<».  Thus, 

x(t|t)  -*■  xk(t  |t)  (3-23) 


in  probability  as  t-*-«®. 

The  other  assumption  which  must  be  relaxed  is  that  the  parameter  set  is 
fixed.  There  are  many  situations  in  which  the  size  of  the  parameter  set  is 
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always  finite  but  increases  with  time,  i.e.,  the  number  of  plausible  hypoth¬ 
eses  can  grow  steadily  with  time.  The  simplest  such  situation  is  when  d«d(t) 
and  d(t)  takes  on  values  from  the  usual  finite,  discrete  set,  . ,Djj} ;  the 

situation  is  similar,  if  more  complex,  when  tracking  multiple  targets.  In 
this  simplest  situation,  there  are  N  hypotheses  at  time  t*l ,  each  correspond¬ 
ing  to  one  possible  value  for  d(l).  At  t*2,  there  are  N2  hypotheses,  each 
corresponding  to  a  plausible  pair  of  values  for  d(l)  and  d(2),  and  so  forth  as 
indicated  in  Fig.  3-2.  Each  node  in  that  figure  represents  a  hypothesis  and 
associated  with  each  hypothesis  is  a  state  estimate  Xi(t|t)  which  is  condi¬ 
tioned  on  the  hypothesis  being  correct  and  the  posterior  probability  of  the 
hypothesis'  correctness,  Pi(t). 

The  conditional  estimate  x^(t|t)  is  the  standard  Kalman  filter  estimate 
for  the  model  of  Eqs.  3-1  through  3-11  with  d(t)=di_[i/N]N>  £f(t)  is  the 
corresponding  filter  residual 

A 

Jji(t)  «  £(t)  -  C(t;d(t))xji/N](t|t-1)  ,  ( 3—24 ) 

Sj[(t)  is  the  residual  covariance  matrix, 

6i(t)  -  jS^t)!1/2  ,  (3-25) 


and 


Pi(t) 


Bj  exp  -  -  rj(t)TSj(t)-1ri(t) 


1 

2 


Nc 

Z  Bj(t)exp  |  -  -  rj(t)TSj(t)-1ri(t)  |  P[j/N](t~l) 


) 


P[i/N](t-1) 


(3-26) 
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d(l)=Dj ,d(2)=Dj 
d(l)=D1,d(2)=D2 


d(l)=DN.d(2)=D1,d(3)=D4,d(4)=D1 


d(t)e{Dlf...,DN}  t=l,2, . . . 

11-069? 


Figure  3-2.  Illustrating  Growing  Trees  of  Hypotheses. 
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The  indices  are  "coded"  so  that  the  first  (i«l)  hypothesis  is  that  d(T)*I)j  for 
Kt,  the  second  (i=2)  is  that  d(T)=D1  for  r<t-l  and  d(t)=JD2,  and  the  last 
(i*Nt)  is  that  <1(t)«Dj;  for  T<t.  The  notation  [i/N]  denotes  the  largest 
integer  n  strictly  less  than  i/N.  For  example,  [ l/2]=[ 1 ]«0,  [4/3]=[2)=l,  etc. 
If  the  dimension  of  the  state  vector  is  not  a  function  of  d,  then  it  is  pos¬ 
sible  to  form  an  unconditional  state  estimate  as  in  Eq.  3-13.  Otherwise,  a 
maximum  a  posteriori  criterion  must  be  used  to  select  one  of  the  conditional 
state  estimates  as  the  state  estimate. 

3.2  OPTIMAL  TRACKING  ALGORITHMS 

As  with  our  model  development  effort,  our  algorithm  development  effort 
has  proceeded  by  considering  three  cases  of  increasing  complexity.  The  cases 
we  have  considered  include: 

1.  one  target  (known)  with  measurement  errors,  missed  detec¬ 
tions,  and  false  detections,  but  no  maneuvers,  etc.; 

2.  one  target  (known)  with  measurement  errors,  missed  detec¬ 
tions,  false  detections,  "mortality,"  (appearance  and 
disappearance)  but  no  maneuvers,  etc.;  and 

3.  an  unknown  number  of  targets  with  measurement  errors, 
missed  detections,  false  detections,  “mortality,"  but  no 
maneuvers,  etc. 

Case  1  is  the  simplest  combination  of  phenomena  which  requires  an  MMAE 
algorithm.  The  only  hypotheses  concern  which  detection  on  a  given  scan  is 
caused  by  one  target,  if  any.  This  is  the  simplest  sort  of  data  association 
hypothesis.  In  Case  2  is  is  necessary  to  hypothesize  whether  the  target  has 
been  "born"  and  if  so,  whether  it  has  also  "died."  These  hypotheses  are,  of 
course,  in  additic-  to  the  simple  hypotheses  about  which  detection,  if  any,  is 
caused  by  the  target.  Case  3  provides  far  more  complexity  in  the  data  asso¬ 
ciation  hypotheses,  as  there  are  many  possible  pairings  of  targets  and 


40 


ALPHATECH,  INC. _ 

detections.  Hypotheses  concerning  the  number  of  targets  born  by  a  given  time 
(as  well  as  which,  if  any,  have  died)  are  also  essential.  Each  case  will  be 
discussed  in  varying  detail  in  the  succeeding  subsections. 

3.2.1  Case  1 

In  this  case  there  is  a  single  target  (and  it  is  known  there  is  just  one) 
which  is  ''immortal"  in  the  sense  that  its  "birth"  and  "death"  are  not  modeled. 
The  unknown  parameters  are  the  single  target's  detection  state,  which  varies 
with  time  and  which  takes  on  one  of  L(t)+1  values  at  time  t.  For  this  reason, 
a  generalized  MMAE  algorithm  of  the  form  described  above  is  needed.  The  algo¬ 
rithm  creates  an  ever  increasing  number  of  hypotheses,  one  corresponding  to 
each  possible  sequence  of  values  that  can  be  taken  on  by  the  sequence  of 
detection  states,  i.e.,  at  time  t  there  is  one  for  each  possible  sequence 
{d(x)}t  . 

All  of  the  information  associated  with  the  i-th  hypothesis  at  t  can  be 
gathered  into  a  global  hypothesis  structure,  denoted  by  Hc(i,t).  It  has  four 
components : 


Hc(i,t)  =  (Dg( i »t ) ,x(t | i ,t ) ,E(t | i , t ) ,PG(i »t ) )  (3-27) 

The  first  component  represents  the  hypothesized  detection  state  sequence.  As 
follows : 


DG(i,t) 


(dG( r|i,t)} 

t— « 


(3-28) 


and 
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d(;(T | i ,t)  -  [bG(T|i,t),  iG(i|i,t)]T  for  x<t 


(3-29) 


Effectively,  dG(T|i,t)  is  the  value  of  d(-r)  assumed  by  the  i-th  hypothesis 
existing  at  t;  bG(T|i,t)  and  iG(x|i,t)  are  similarly  related  to  b(x)  and  i(x). 
The  i-th  hypothesis  at  t  is  correct  if  and  only  if 


dc(x|i,t)  =  d(x)  for  all  x<t 


(3-30) 


For  convenience,  we  denote  this  event  by  ©G(i,t). 

The  second  and  third  components  of  HG(i,t)  represent  the  optimal  estimate 
of  the  single  target's  dynamic  state  x(t),  as  follows: 


x(  1 1  i  ,t)  »  E[x(t)|0G(i,t)&ZG(t); 


(3-31) 


and 


I(t|i,t)  -  Cov [ x( t ) | 0G( i,t)&ZG(t)J 


(3-32) 


where  E[ • | • ]  and  Cov[*|*]  are  respectively  the  conditional  mean  and  condi¬ 
tional  covariance  operators,  and  where 


ZG(t)  -  {Zs( t) }  for  all  t 

T*-® 


(3-33) 


and 


L(t) 

Zs( t)  »  {z(j,x)}^  ^  for  all  x<t 


(3-34) 
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The  Zs(t)  sequences  include  all  of  the  detections  for  each  scan  t  while  the 
ZG(t)  sequences  include  all  of  the  detections  in  all  of  the  scans  up  through 
t.  Since  x(t)  is  Gaussian  distributed  for  all  t  (see  Eq.  2-1),  x(t|i,t)  and 
I(t|i,t)  together  give  the  conditional  distribution  of  x(t).  Note  that  the 
definitions  of  x(t|i,t)  and  I(t|i,t)  can  be  extended  to  making  x(i|i,t)  and 
r ( t  |  i  , t )  the  equivalent  conditional  mean  and  covariance  of  x(x)  for  T*t. 

The  fourth  component  of  H^Ci.t)  represents  the  probability  that  the 
hypothesis  is  supported  by  the  measurement  data  as  follows: 

PG(i,t)  =  Prob[ 0G( i ,t ) | ZG( t ) ]  .  (3-35) 

We  denote  the  number  of  hypotheses  at  t  by  M(;(t). 

The  hypothesis  structures  are  the  data  structures  manipulated  by  the  op¬ 
timal  algorithm.  It  is  necessary  to  define  how  the  components  are  calculated. 
Before  doing  so,  it  is  convenient  to  define  one  more  data  structure,  the  scan 
hypothesis  structure,  denoted  by  Hs(j,t).  It  has  one  component: 

Hs(j,t)  '  (ds(t|j,t))  for  all  t  (3-36) 


where 


=  IbsCc  |j  »t>,  isCt  |j  ,t)JT  for  all  t  .  (3-37) 


Effectively,  ds(t|j,t)  is  the  value  of  d(t)  assumed  by  the  j-th  scan  hypoth¬ 
esis  at  t;  bs(t|j,t)  and  is(t|j,t)  are  similarly  related  to  b(t)  and  i(t). 
The  j-th  scan  hypothesis  is  correct  if  and  only  if 
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ds(t|j,t)  -  d(t)  .  (3-38; 

For  convenience,  we  denote  this  event  by  0s(j,t).  We  denote  the  number  of 
scan  hypotheses  at  t  by  Ms(t).  In  this  case,  the  only  plausible  scan  hypoth¬ 
eses  are  that  the  single  target  caused  one  of  the  L( t )  detections  or  that  it 
was  not  detected,  so 

Ms(t)  =  L(t)  +  1  for  all  t  .  (3-39) 

Each  global  hypothesis  structure  is  constructed  from  a  sequence  of  con¬ 
secutive  scan  hypothesis  structures.  Thus,  for  each  global  hypothesis  struc¬ 
ture  at  t-1  and  scan  hypothesis  structure  at  t,  one  can  construct  a  new  global 
hypothesis  structure  for  t.  Symbollically , 

HG(i,t-l)*Hs(j,t)  =>  Hc(k,t)  (3-40) 

for  every  i  and  j  pair.  It  is  appropriate  to  note  at  this  point  that  the  in¬ 
dices  of  global  and  scan  hypothesis  structures  are  arbitrary.  One  could 
cleverly  encode  the  value  of  DG(i,t)  in  the  index  i  and  of  ds( j ,t)  in  the 
index  j,  making  the  hypothesis  structures  "content  addressable."  But  the 
choice  of  whether  to  do  so  or  not  is  basically  an  implementation  issue  which 
we  wish  to  avoid.  In  defining  the  optimal  tracking  algorithm,  we  will  simply 
state  it  generally  in  terms  of  arbitrary  indices  i,  j,  and  k  as  in  Eq.  3-40 
without  stating  how  k  should  be  chosen  for  a  particular  i  and  j  pair. 

The  new  hypothesized  detection  state  sequence  DG(k,t)  is  generated  from 
DG(i,t-l)  and  ds(j ,t)  by  letting 

dG(t|k,t)  «  dG(r|i,t-l)  for  x<t-l  (3-41a) 
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and 


dG(t  |k,t)  -  ds(  t  |j  ,t) 


(3-Alb) 


Assuming  that  bs(t|j,t)=D,  the  dynamic  state  estimate  is  updated  by  letting 


x(t |k,t)  =  x( t ) i , t-1 )  +  K( 1 1 i ,t-l ) [ z(  is(t | j ,t ) ,t )-Cx( t | i , t-1 ) ]  (3-Alc) 

and 

E(t |k,t)  =  H-K(t|i,t-l)CTjE(t|i,t-l)  .  (3-Ald) 

If  bs(t|j,t)=D,  there  is  no  measurement  to  incorporate,  ad 

x(t  jk,t)  =  x(t  j i ,t-l )  .  (3-Ale) 

Here 

x(t|i,t-l)  =  Ax( t-1 | i , t-1 )  ,  ( 3-A2a) 


E(t | i  ,t-l )  =  AI(t-l | i , t - 1 ) At  +  W  , 


( 3— A  2b) 


K( t  | i ,t-l )  -  Z(t|i,t-l)CT[v+CE(t|i,t-l)CTj-1  . 


(3-A2c) 


The  probability  is  updated  by  letting 


PG(k,t)  -  PG(i,t-l)(l-PD)Q(t)  if  bs(t|j,t)  =  D  (3-A3a) 


PG(i,t-l)PD[G(z(is(t|j,t),t)-x(t|i,t-l),V+CE(t |i,t-l)CT)/pFA)Q(t) 

a 

if  bs ( t )  j  ,  t )  -  D 
A5 


( 3-A3b) 
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where  G(*,»)  is  the  Gaussian  probability  density  function,  i.e., 

-yTy-ly/2 

C(y,Y)  -  (2w)-K/2  j  Y  j-l/2e  ,  (3-44) 

and 

"*FA  L(t) 

Q(t)  =  [e  [pFA  /L( t ) ! ] j/Prob[Zs(t) | ZG( t-1 ) ]  .  (3-45) 

This  algorithm  is  carried  out  for  every  pair  of  i  and  j  such  that  l<i<M^(t-l) 
and  l<j<Ms(t).  It  follows  that 

Mc(t)  =  Mg( t-1 )MS( t )  =  Mc(t-l) [L(t )+l J  for  all  t  ,  (3-46) 

The  practical  infeasibility  of  the  algorithm  is  here  made  clear;  the  number  of 
global  hypotheses  grows  geometrically  with  t. 

Note  that  it  would  be  possible  to  form  an  unconditional  state  estimate 

iicCt )  . 

x(t|t)  =  Z  x(t  |k,t)PG(k,t)  (3-47) 

k=l 

in  this  case  because  the  dimension  of  x  is  always  4.  But  we  assume  for  con¬ 
sistency  with  more  complex  cases  that  a  maximum  a  posteriori  criterion  is 
used.  That  is,  for  each  t 

A  A 

x  ( 1 1 1 )  =  x( 1 1 k* , t )  (3-48) 

where  k*  is  such  that 

PG(k*,t)  >  PG(k,t)  for  all  k  .  (3-49) 


F 
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The  equations  generating  the  new  hypothetical  detection  state  sequence 
are  obviously  appropriate.  The  equations  generating  the  new  dynamic  state  es¬ 
timate  are  nearly  as  obvious;  they  are  the  Kalman  filter  equations  specialized 
to  the  particular  pairing  of  hypotheses.  This  is  the  optimum  way  to  update  a 
dynamic  state  estimate  for  a  Gaussian  distributed  dynamic  state,  given  a  mea¬ 
surement  based  on  that  dynamic  state.  The  equations  generating  the  probabil¬ 
ity  of  the  new  global  hypothesis  are  less  obvious.  Their  derivation  is 
sketched  in  the  next  few  paragraphs. 

By  definition, 

PG(k,t)  =  Prob[ 0G(k , t ) j  Zg( t ) ]  .  (3-50) 

Observing  that  the  event  0G(k,t)  occurs  if  and  only  if  the  events  0G(i,t-l) 
and  0s(j,t)  occur  and  that 

ZG(t)  =  ZG(t-l)U{Zs(t)}  ,  (3-51) 


one  can  demonstrate  that 

PG(k,t)  =  Prob[Zs(t) | 0S(  j , t )&0G( i ,t-l )&ZG( t-1 ) ] 

x  Prob[Os(j , t ) | ©G( i , t— 1 ) &ZG( t  —  1 ) J  (3-52) 

X  PG(i,t-l)/Prob[Zs(t) I ZG( t— 1 ) ] 

The  first  factor  on  the  right-hand  side  of  the  above  equation  can  be  further 


expanded  as  follows 
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Prob [ Zs( t ) | 0S( j ,t)&eG(i,t-l)&ZG(t-l)] 


L(t) 

-  n  Prob [ z( 4  ,t  ) | L ( t ) & ©s ( j , t )&0G(k , t-1 )&ZG( t-1 ) ] 
4=1 


x  Probj  L( t ) | 0S( j ,t )&0c(i ,t-l )&ZG( t-1 ) ]  (3-53) 

because  the  measured  locations  of  all  detections  are  statistically  indepen¬ 
dent.  Further  simplification  of  Eqs.  3-52  and  3-53  can  be  accomplished  by 
noting  similar  independence,  leading  to  the  observation  that 

L(t) 

PG(k,t)  =  n  Prob[zU,t)|0s(j,t)&0G(i,t-l)&ZG(t-l)] 

4=1 


x  Prob [ L( t ) | 0S( j ,t ) ] Prob [ 0S( j , t) J PG( i , t-1 )/ Prob [ Zs( t ) | ZG( t-1) ] 

(3-54) 

Further  simplification  requires  assumptions  about  Hs(j,t). 

Assume  that  j  is  such  that  ds(j,t)=D.  This  implies  all  detections  are 
false  alarms.  The  measurements  of  all  of  the  detections  are  uniformly  dis¬ 
tributed  and  LpA(t)  equals  L(t).  Thus 

L(t)  -*FA  L(t) 

PG(k,t)  =  [  n  |=|_1][e  XFA  /L(t ) ! J [ 1-PD] PG(i , t-1 )/Prob[ Zs( t ) | ZG( t-1 ) J  . 
4=1 

(3-55) 

A  little  manipulation  yields  Eq.  3-43a. 

Now  instead  assume  that  j  is  such  that  ds(t|j,t)=D.  This  implies  that 
all  detections  but  the  is(t|j,t)-th  are  false  alarms.  The  measurements  of  all 
the  false  detections  are  uniformly  distributed  and  LpA(t)  is  one  less  than 
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L(t).  The  measurement  z(  i(t | j ,t)  ,t)  is  Gaussian  distributed  if  0s(j,t).  From 
Eq.  2-12a  it  is  clear  that  the  mean  and  covariance  of  the  measurement,  condi¬ 
tioned  as  above  are 

Cx( t | i , t-1 )  (3-56a) 

and 

V  +  CZ(t |i,t-l)CT  ( 3-56b)  ' 


respectively.  It  follows  that 

Prob(z( is(t | j ,t),t) |0s(j ,t)&0G(i ,t-l )&ZGC  t-1 ) J 

=  G(z(xs(t|j,t),t)-Cx(t|i,t-l),V+CE(t|i,t-l)CT)  (3-57) 


and  sc 


L(t) 

PG(k,t)  =  [  n  |~|-1]G(z(is(t|j,t) ,t)-Cx( t | i , t— 1 ) ,V+CI(t | i , t — 1 ) CT) 

1=  1 

i*is(t|j  ,t) 


~XFA  L(t)-1 

^  [e  [XFA  ]/[L(t)-l]!]( PD/ L(t) )PG(i,t-l)/Prob[Zs(t) |ZG(t-l) ] 

(3-58) 


A  little  manipulation  yields  Eq.  3-43b. 


3.2.2  Case  2 

In  this  case,  there  is  again  a  simple  target  (and  it  is  again  known  there 
i 8  just  one),  but  it  is  "mortal"  in  that  its  "birth"  and  "death"  are  modeled. 
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The  unknown  parameters  are  again  the  target's  detection  state,  which  varies 
with  time  and  which  takes  on  one  of  L(t)+3  values  at  time  t.  Moreover,  the 
dynamic  state  of  the  target  is  effectively  of  dimension  0  (i.e.,  it  is  un¬ 
defined)  if  d(t)=B  or  X.  For  both  of  these  reasons,  a  generalized  MMAE  algo¬ 
rithm  of  the  form  described  in  subsection  3.1.2  is  needed.  Again,  the  algo¬ 
rithm  creates  an  ever  increasing  number  of  hypotheses,  one  corresponding  to 
each  possible  sequence  of  values  that  can  be  taken  on  by  the  sequence  of 
detection  states,  i.e.,  at  time  t,  there  is  one  for  each  possible  sequence 
{d(  t)  }  t=q  . 

The  modifications  to  the  model  going  from  Case  1  to  Case  2  are  not  exten¬ 
sive;  the  modifications  of  the  optimal  tracking  algorithm  are  also  not  exten¬ 
sive.  The  global  and  scan  hypothesis  structures  are  defined  exactly  as  in 
Case  1.  Only  the  range  of  values  of  bG(i,t)  and  bs(i,t)  differ.  We  can  again 
symbolize  the  generation  of  new  global  hypothesis  structures  for  old  global 
hypothesis  structures  and  current  scan  hypotheses  by 

HG(i,t-l)*Hs(j,t)  =>  HcCk.t)  (3-59) 

but  with  two  restrictions.  The  first  restriction  is  that  this  generative  pro¬ 
cedure  applies  only  if  t>0.  For  t=0,  we  have  but  one  global  hypothesis  struc¬ 
ture,  Hg(1,0),  such  that 

bG( 1 ,0)  -  B  (3-60) 

and 

PG(1,0)  -  1  .  (3-61) 
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All  of  the  other  components  of  this  global  hypothesis  structure  are 
undefined.  The  second  restriction  is  that  Eq.  3-59  does  not  apply  for  every  i 
and  j  pair.  It  applies  only  if  the  event  state  ds(j,t)  is  reachable  from 
dG(i,t-l)  in  one  step,  i.e.,  if  either 

bs  ( c  I J  » t )  =  B  and  bG( t-1 | i ,t-l )  =  B  , 

bs( 1 1  j ,t)  =  D  and  bG( t-2 j i  , t-1 )  *  X  , 

bs  C c  I J  » c )  “  D  and  bG( t-1 | i ,t-l )  *  D  or  D  ,  or 

Mtb.t)  “  X  and  bG(  t-1 1  i ,  t-1 )  =  D  or  X  .  (3-62) 

Given  the  aforementioned  restrictions,  we  can  describe  the  generative 
procedure.  The  new  hypothetical  event  state  sequence  D^k.t)  is  generated 

A  A 

from  DG(i,t-l),  and  ds(j,t)  by  letting 

dG(t|k,t)  *=  dG(  x|  i  ,t-l )  for  0<Kt-l  (3-63a) 

and 


! 


dG(t|k,t)  -  ds(t  |  j  ,t)  ,  (3-63b) 

exactly  as  in  Case  1. 

Updating  the  dynamic  state  is  a  bit  more  complex  in  this  case,  the  for¬ 
mulae  depending  on  the  hypothesized  event  state  in  a  manner  paralleling  Eq. 

A  A 

2-21.  If  bs(t|j,t)«D  and  bc<t-l |i,t-l)-D  or  D, 
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x(c|k,t)  -  x(t|i,t-l)  +  K(t|i,t-l)(z(ts(t|j ,c) ,t )-x( t |i,t-l) J  (3-64.) 

and 


X(t|k,t)  -  lI-K(t|i,t-l)cT]r(t|i,t-l) 


while  if  bs(t|j,t)«D,  then 


x(t|k,t)  =  x(t|i,t-l) 


and 


r(c  jk.t)  =  x(t|i,t-i) 


where 


x(t|i,t-l)  =  Ax( t-1 | i , t-1)  , 


E( 1 1  i  ,  t-1 )  =  AZ(t-l ji,t-l)AT  +  W  , 


and 


But  if  bs(t|j,t)-D  and  bG( t-1 | i , t-1 )=B ,  then 


x(t|k,t)  -  KqZ( ls(t | j ,t) ,t) 


( 3-64b) 


( 3-65a) 


(3-65b) 


(3-66a) 

(3-66b) 


K(t|i,t-1)  =  E(t|i,t-l)cT|V+CZ(tji,t-l)cT]-l  .  (3-66c) 


(3-67a) 


and 


Z(t|k,t)  -  [I-KqCTJZc 


( 3—67  b ) 
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where 

Ko  *  E0Ct[V+CE0Ct]-1  .  (3-68) 

For  all  other  circumstances,  x(t|k,t)  and  Z(t|k,t)  are  undefined. 

The  updating  of  the  probability  is  similarly  dependent  on  the  values  of 
bs(t  |  j  ,t)  and  bG(t-l | i ,t-l ) . 

PG(k,t)  -  PG( i ,t— 1) [ 1— PB]Q(t)  if  bs(t|j,t)=B  and  bG(t-l | i ,t-l )-B 

PG(i,t-l)PB[G(z( is(t|j,t),t),V+CZ0CT)/pFA]Q(t) 

A  A 

if  bs(t|j,t)=D  and  bG(t-l | i ,t-l)=B 
PG(i,t-l)PD[G(z(is(t|j,t),t)-Cx(t|i,t-l),V+CZ<t|i,t-l)CT)/pFA]Q(t) 


if 

bs(t|j,t)-D 

and 

bG( t-1 | i ,t-l )=d  or  D 

PG(i,t-l)[l-PD-Px]q(t) 

if 

A 

bs(t|j,t)=D 

and 

A 

bG( t-1 | i , t-1 )=D 

PG(i,t-l)PXQ(t) 

if 

bs(t|j,t)-X 

and 

bG(t-l | i ,t-l )=D 

PG(i ,t-l ) [ 1— PD] Q( t ) 

if 

bsC  1 1 j ,t)-D 

and 

bG( t-1 )i ,t-l )=D  , 

PG(i,t-l)Q(t) 

if 

bs(t|j,t)-X 

and 

bG( t-1 | i , t-1 )=X  . 

-yTy-ly/2 

G(y,Y)  -  C  2  it )  2  |  y  j-l/2e  (3-70) 
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and 


“*FA  L(t) 

Q(t)  =  [e  [pFA  /L(t)! j]/Prob(Zs(t) |ZG(t-l)]  (3-71) 


as  before. 

In  this  case,  the  plausible  scan  hypotheses  at  t  are  that  the  single  tar¬ 
get  caused  one  of  the  L(t)  detections,  that  it  has  not  been  "born"  yet,  that 
it  is  already  "dead"  or  that  it  simply  is  undetected  but  not  yet  "dead."  So, 

Ms ( t )  -  L(t)  +  3  for  t>0  .  (3-72) 

Given  the  observation  that  not  every  combination  of  an  old  global  hypothesis 
and  a  current  scan  hypothesis  generates  a  new  global  hypothesis,  the  most  that 
can  be  said  is  that 


flG(t)  <  Mg( t-1 )  [ L(t )+3]  for  t>0  ,  (3-73) 

meaning  that  the  growth  of  the  number  of  global  hypotheses  can  be  geometric. 

The  derivation  of  the  probability  update  formulae  is  quite  similar  to 
that  in  Case  1.  One  difference  appears  in  the  equivalent  of  Eq.  3-34.  In 
this  case, 

L(t) 

PG(k,t)  -  n  Prob[z(£,t)|os(j,t)&0G(i,t-l)&ZG(t-l)] 

t-1 

x  Prob[L(t)|0s(j ,t)]Prob[0s(j , t) | 0G( i ,t-l) ] PG(i ,t-l )/Prob[ Zs( t ) | ZG( t-1 ) ]  . 

(3-74) 
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The  third  factor  differs  from  Eq.  3-54;  the  probability  of  a  hypothesized 
event  state  at  t  depends  on  the  hypothesized  event  6tate  at  t-1.  Another  dif¬ 
ference  appears  in  the  equivalent  of  Eq.  3-57  when  bs(t|j,t)*D  and 
b(j(  t-1  |  i  ,  t  )“B.  The  measurement  z(  is(  t  j  j  , t ) ,  t )  is  Gaussian  distributed  if 
Hs(j,t)  is  true,  but  its  mean  and  covariance  are,  respectively 

CE(x0]  -  0  ( 3-75a) 

and 

V  +  C  Cov  [  x0 ]  CT  ■=  V  +  CI0CT  (3-75b) 

respectively.  It  follows  that 


Prob[z(is(t|j,t),t)|0s(j  ,t)&0G(i  ,t-l  )&ZG(  t-1 )  )  =  G(z(is(t|j,t),t),V+C5;0CT)  . 

(3-76) 

3.2.3  Case  3 

In  this  case,  there  are  an  unknown  number  of  targets,  each  of  which  can 
be  "born"  or  "die"  at  any  time.  The  unknown  parameters  include  the  number  of 
targets  having  been  born  and  their  detection  states,  each  of  which  can  take  on 
L(t)+2  values.  The  overall  dynamic  state  vector  is  composed  of  the  four 
dimensional  dynamic  state  vectors  of  the  targets  which  are  currently  "alive," 
making  the  dimension  of  the  overall  dynamic  state  vector  a  function  of  the  un¬ 
known  parameters.  For  both  of  these  reasons,  a  generalized  MMAE  algorithm  of 
the  form  described  in  subsection  3.1.2  is  needed.  Again,  the  algorithm 
creates  an  ever  increasing  number  of  hypotheses,  one  corresponding  to  each 
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possible  sequence  of  values  which  can  be  taken  on  by  the  sequence  of  composite 

discrete  states,  i.e.,  at  time  t,  there  is  one  for  each  possible  sequence 
t 

{d(t)}T,o  where 

d(i>  -  [N(T),d'(l,t)T,d'(2,t)T . d'(N(t),t)T]T  (3-77) 

as  in  subsection  2.1.3. 

The  global  and  scan  hypotheses  are  defined  exactly  as  in  Case  1,  but  the 
variables  d(t)  and  x(t)  now  have  new  components.  Hypotheses  about  and 
estimates  of  these  new  components  are  denoted  as  usual.  Again  we  can  symbol¬ 
ize  the  generation  of  new  global  hypothesis  structures  from  old  global  hypoth¬ 
esis  structures  and  current  scan  hypotheses  by 

HG(i,t-l)*Hs(j,t)  =>  HcCk.t)  (3-78) 

but  with  restriction,  as  we  did  in  Case  2.  The  restrictions  are  fewer,  but 
slightly  more  complex.  For  t“0 ,  we  again  have  but  one  global  hypothesis 
structure,  Hc(l,0).  Its  components  are 

NG(1,0)  •=  0  (3-79) 

and 

PG(1,0)  -  1  .  (3-80) 

The  second  restriction  is  that  Eq.  3-77  applies  only  if  the  total  event  state 
ds(j,t)  is  reachable  from  dG(i,t-l)  in  one  step,  i.e.,  only  if 

Ns(j,t)  >  NG(i,t-l)  (3-81) 
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and  if  for  every  l  either 


bsU,t|j,t)  -  B  and  bG(  l ,  t-1  |  i  ,  t-1 )  -  B 


(3-82a) 


bs(£,t | j ,t)  -  D  and  bG( i »t-l | i ,t-l )  *  X  , 


( 3-82b) 


bs(i,t  |  j  ,t)  -  D  and  bG( l ,t-l | i ,t-l )  -  D  or  D 


(3-82c) 


bsU,t|j,t)  *  x  and  bG( 4,t-l | i , t-1 )  =  D  or  X 


( 3-82d ) 


Some  additional  notation  is  needed  before  the  generation  procedure  can  be 
derived.  We  partition  the  set  of  target  indices  as  follows: 

a  A  A  ___ 

H'DB(k,t)  «=  U|bsU,t|j,t)=D&bGU,t-l|j,t-l)=B}  (3-83a) 

fDD^.b)  =  |bsU,t  |  j  ,t)»D&bGU,t-l  |j  ,t-l)“D}  (3-83b) 

^DD(k.t)  -  U|bsU,t|j,t)«D&bGU,t-l|j,t-l)«D}  (3-83c) 

^XD(^.t)  “  {£■  |  bs(  4 1 1 1  j  » t  )“X&bG(  i,t-l  |i,t-l)»D}  (3-83d) 

A  A  A 

fDD(k,t)  “  {t|bs(i,t|j,t)-D4bc(i.t-l|j,t-l)-D}  (3-83e) 


“  U|bs(*.t|jft)-D&bG(*,t-l|j,t-l)-D}  (3-83f ) 
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txx(k.t)  =  U|bsU,t)j,t)«X&bGU,t-l  jj,t-l)=X}  .  (3-83g) 

The  cardinality  of  each  of  these  sets  is  denoted  by  using  the  corresponding 
lower  case  symbol,  e.g., 


itoB(k.,t)  *=  Card!  fDB(k,t)  ]  (3-84) 

The  new  hypothesized  total  event  state  sequence  DgCk.t)  is  generated  from 
Dgdjt-l)  and  ds(j,t)  by  letting 

dG(  t  |k,t )  =  dG(i|i,t-l)  for  0<t<t-l  (3-85a) 


and 


dG(t |k,t)  *  ds(t  |  j  ,t) 


(3-85b) 


exactly  as  in  Case  1  but  for  the  restrictions  on  i  and  j. 

Updating  the  dynamic  state  in  this  case  resembles  the  process  in  Case  2: 


x* ( A , 1 1 k, t )  =  x ' ( £ , t  j t-l)+K( £,t|i,t-l) I z( is(£,t | j ,t) ,t)-Cx' ( £,t J i , t  — 1 ) ] 

(3-86 a) 


and 


r(i,t|k,t)  =  [ I-K( £,t|i,t-l)CT]E'(£,t|i,t-l)  (3-86b) 

for  tc’rDD(k,t)UyDD(k,t)  and 


x'(£,t|k,t)  *  x'(£,t | i , t— 1 ) 


(3-87a) 
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and 


EU,t|k,t)  =  (3-87b) 

for  £ef[)[)(k ,t )U'll[)D(k )  where 

x'(£,t | i , t — 1 )  =  Ax'U.t-l |i,t-l)  ( 3-88a ) 

Z'(£,t|i,t-1)  =  AZ’U.t-l  |i,t-l)AT+W  (3-88b) 

and 


K(A,t  | i , t — 1 )  =  r'U.tli.t-ncTlV+CI'Uli.t-ncT)-1  ,  (3-88c) 


plus 


x'U,t|k,t)  =  KqZ(  is(Jt,t  |  j  ,t),t) 


(3-89a) 


I’U,t|k,t)  =  II-KoCt]e'  (3-89b) 

o 

for  £e¥[)B(k,t)  where 

kg  -  r'cT[v+cr'cT)-i  .  (3-90) 

O  O 

Under  all  other  circumstances,  x'(£,t|k,t)  and  E'(i,t|k,t)  are  undefined. 

Updating  the  probability  more  extensively  involves  the  partitions  of  the 
target  indices: 
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PG(k,t)  -  PG(i,t-l)  n  XB[G(z(is(i,t!j,t),V+CE’CT)/pFA]  x 

n  PD[G(*(i8(A,t|j.t) ,t )-Cx( £)  ,  1 1  i  ,  t-1 ) 

teVk*t)uWk*t) 

V  +  CZ'U.t  |i,t-l)C'T)/pFA]  X 

^DD(k>t)  *XD(k)t)  Wk*t} 

(1-PD-PX)  PX  (1-Pd)  Q(t)  (3-91) 


where 


_yTY-ly/2 

G(y,Y)  =  (2tt)_k:/ 2  |  Y  |~l/2e 


(3-92) 


and 


"*FA  L(t) 

Q(t)  =  e  [e  [ pFA  / L( t ) ! ] ]/Prob[Zs(t)/ZG(t-l) ]  .  (3-93) 

Note  that  Q(t)  does  differ  from  the  previous  case,  but  still  is  independent  of 
the  hypothesis  indices. 

In  this  case,  the  plausible  scan  hypotheses  at  t  are  quite  large.  As 
many  as  L(i)  new  targets  or  as  few  as  0  could  have  been  born  on  each  scan  T<t, 
so  one  knows  only  that 


t 

0  <  N(t)  <  Z  L(t)  .  (3-94) 

t-1 


60 


ALPHATECH,  INC 


In  terms  of  the  Ns(t)  components  of  the  scan  hypothesis  structure  alone,  there 
can  be 


t 

Z  L(t) 

i=  1 


+  1 


(3-95) 


distinct  scan  hypotheses.  The  various  possiole  values  of  ds(t|j,t)  are 
naturally  larger  but  can  be  bounded  loosely  as  follows: 


Ms(t)  <  [L(t)+2j 


t 

I  L(t) 

T—  1 


+  1 


for  t>0 


(3-96) 


Because  the  observation  that  not  every  combination  of  an  old  global  hypothesis 
and  a  current  scan  hypothesis  generates  a  new  global  hypothesis  is  especially 
true  in  this  case,  it  can  be  fairly  said  that 


MgCO  <  Mc(t-l)[L(t)+2] 


t 

Z  L(t) 

T=1 


+  1 


for  t>0 


(3-97) 


is  a  very  weak  bound.  This  bound  does  reflect  a  potentially  faster  that  geo¬ 
metric  growth  of  the  number  of  hypotheses  in  this  case. 

The  derivation  of  the  probability  update  formulae  is  naturally  more  com¬ 
plex  in  this  case,  but  does  parallel  that  of  Cases  1  and  2.  Equations  3-52 
and  3-53  are  still  valid,  for  example.  The  latter  equation's  terms  can  be 
further  expanded  as  follows: 
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L(t) 

n  Prob[z(A,t) i L( t )&0S( j ,t)&0G(i,t-l)&ZG(t-I) j 
4=1 


-  -  n  *  Prob[ z( is( A , t | j , t ) , t ) | 0S( j , t )&©G( i , t-1 )&ZG( t-1 ) ] 

AeT  _( k  ,  t  )l'¥  (k,t)Uf  _(k,t) 

DB  DD  DD 


Prob[z(A,t)]  (3-98) 

Ae4>(k,t) 


and 


Prob[L(t)|0s(j,t)&0G(i,t-l)&ZG(t-l)]  =  Prob[  <f(k  ,t )  j  (3-99) 


where 


*(k,t)  =  {  1  , . . .  ,L(  t ) }  -  {ls(  A,t  |  j  ,t)  |  AE^B(k»t)u^DD(k;»t)u,toD(k>t)  (3-100) 


is  the  set  of  indices  of  the  detections  implicitly  hypothesized  to  be  false, 
and 


<t>(k  ,t )  =  Card  [  $(k  ,t )  j 


(3-101 ) 


Reasoning  about  the  conditional  means  and  variances  of  the  measurements  re¬ 
sulting  in  target  detections  as  in  Cases  1  and  2  leads  to  the  conclusion  that 
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Prob[Zs(t) |0s(j ,t)&eG(l,t-l)&ZG(t-l)]  - 


*  n  G(z(is(t,t|j,t),t),V+CZ,cT)  x 

fet  _(k,t) 

DB 


11  -  G(z(  is(£,t |j,t),t)-Cx'(£,t |i,t-l), 

lev  (k.t)U'F  _(k,t) 

DD  DD 


V+CST(£,t|i,t-l)CT)x 


-<Kt)  -ApA  $(t)  * 

e  AFA  /  <#»(  t ) !  .  (3-102) 


The  observation  that 


4>(  t )  =  L(t)  -  ~  (3-103) 

is  used  later  to  eliminate  $(t)  from  the  above  equation. 

Equation  3-102  gives  a  detailed  expansion  of  the  first  term  of  Eq.  3-52 
for  this  case.  The  expansion  of  the  second  term  is  sketched  next.  First,  we 
note  that 

Prob[6s(j ,t) | GG( t , t-1 )&ZG( t-1 ) ]  = 

Prob[dU,t)=dU,t|j,t)  for  *<N(  t)  |  N(t  )-Ns(j  ,t)&0G(i  ,t-l)&ZG(t-l) )  x 

A 

Prob(N(t)»Ns( j ,t) |©G(i , t-1 )&ZG( t-1 ) ]  .  (3-104) 
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The  second  factor  can  be  manipulated  to  show  that 

*  “*B  _ 

Prob[N( t )mNs( j  ,t )  |  0G(i  ,t-l )&ZG( t-1)  ]  =  e  XB  /«Idb(  A.t)  i  (3-105) 

by  observing  that  the  number  of  targets  "born'  in  the  scan  hypothesis  struc¬ 
ture  simply  equals  the  number  already  "born"  in  the  old  global  hypothesis 
structure  plus  those  just  born. 

The  firsc  factor  on  the  right-hand  side  of  Eq.  3-104  can  be  further  ex¬ 
panded  and  simplified  to  yield 

Prob[d(£,t)=d(£,t|j  ,t)  for  £<H(t) )N(t)=Ns(j , t )&0G( i , t-1 )&ZG( t-1 ) ] 

A  AAA 

«  Prob[  iU,t)  =  isU,t  |  j  ,t)  for  1  ^ddOc  > t  JUH'ddU  )  I 

A  A  A 

^DB(  £  )  =  ^DB(^  >  c  c  )  =  '*,DD(k,t)]  * 

Ng(£,t) 

n  Prob [ b( £ , t )=bs( £,t | j , t ) | N(t )=NS ( j ,t)& 

£=  1 

A  A 

N(t-l)  =  NG(i ,t-l )&b( £, t-1 )=bG( A, t-1 1 i , t  —  1 ) ]  (3-106) 

The  number  of  different  ways  ^B( k » t » t)+'/DD^ k c )  indices  can  be  "drawn” 
from  the  L(t)  measurement  indices  is 

L(t )  !  /  [  L(t)-\|q3B(t  )“'l'DD(t)~i|'DD(t )  1  !  (3-107  ) 
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But  the  requirement  that  the  indices  in  'i'DBtk.t)  be  ordered  by  the  magnitudes 
of  the  corresponding  measurements  reduces  that  quantity  by  >(T)B(t)!,  so  the 
first  factor  of  Eq.  3-106  can  be  simplified  as  follows: 


Probf  i(4,t)=is(£,t|j,t)  for  iefDB(k,t)U'{lDD(k,t)U'l'DD(k,t)  | 


*DB<0  -  fDB(k,t)&'}'DD(t)  =  fDD(k,t)6fDD(t)  =  fDD(k,t)] 


-  [L(t)-i|q)B(t)-tDD(t)-t(T)D(t)]!toB(t)!/L(t)!  (3-108) 


The  remaining  factors  of  Eq.  3-106  are  simply  state  transition  probabilities, 
each  repeated  as  many  times  as  a  particular  transition  is  hypothesized  to 
occur.  Thus, 


Prob[  0S(  j  ,  t)  |  0c(i  ,t-l)&Zc(  t-1)  ]  =(  [L(t  )“'l'DD(t )  1  !  toB(t ) !  /L(t ) !  )  * 


4DD(k»t)  ta(k.t) 

Pb  U-Pd-PxJ  px  pd  H-pd)  1 


~^B  >(DB(k>t)  _ 

e  Xb  / 'f'DB^  t ) !  (3-109) 


When  this  equation  and  Eq.  3-102  are  substituted  in  Eq.  3-32,  the  formula  for 
P(;(k,t)  can  be  obtained  by  a  little  rearrangement  of  terms. 
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3.2.4  Case  3  Revisited 

A  remarkable  feature  of  the  formula  in  Eq.  3-91  for  computing  the  proba¬ 
bility  of  the  k-th  hypothesis  at  t  is  that  the  actual  number  of  detections 
appears  only  in  the  term  Q(t),  common  to  all  of  the  hypotheses  at  that  time. 
While  different  hypotheses  may  implicitly  involve  hypothesizing  different  num¬ 
bers  of  false  detections,  the  differing  numbers  do  not  appear  explicitly;  the 
weighting  of  each  Gaussian  density  function  by  the  inverse  of  pfa  effectively 
accounts  for  the  false  detections.  Another  remarkable  feature  is  the  presence 
of  a  term  corresponding  to  each  target  hypothesized  to  have  been  born  at  that 
time.  This  suggests  a  different  organization  of  hypothesis  structures  and  of 
the  optimal  algorithm  which  may  be  more  readily  Implemented  in  a  practical 
manner . 

The  first  step  in  defining  this  alternative  approach  is  to  define  the 
target  hypothesis  structure,  Hx(i,t).  It  has  four  components: 

Ht( i ,t )  =  (d' (i , t ) ,x( t | i ,t ) ,Z( t | i , t ) ,PT( i ,t ) )  .  (3-110) 

T 

The  first  component  represents  a  hypothesized  event  state  sequence  for  some 
target,  index  unspecified,  as  follows: 

*  A  £ 

DT(i,t)  *  {dT(T|i,t)}  for  t>0  (3-111) 

T-0 


and 


d ' ( t  |  i  , t )  -  I bx( T  j i , t ) ,  ix( t | i ,t ) ]T  for  T<t  (3-112) 

T 
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and  where  dx(T|i,t)  is  the  value  of  d'(4,T)  assumed  by  the  i-th  target  hypoth¬ 
esis  for  some  unspecified  4;  t>x(T|i,t)  and  ix(i|i,t)  are  similarly  related  to 
b( 4 , t)  and  i(4,t).  The  i-th  target  hypothesis  at  t  is  correct  if  and  only  if 

dx(i|i,t)  =  d'(4,T)  for  0<T<t  and  some  4  .  (3-113) 

For  convenience,  we  denote  this  event  by  0^(1 >t| 4). 

The  second  and  third  components  of  Hf(i,t)  represent  the  optimal  estimate 
of  the  target's  dynamic  state  x(4,t),  as  follows: 

xT(t|i,t)  =  E[x(4,t)!er(i,t|4)&2G(t)J  (3-114) 

and 

ET(t|i,t)  =  Cov [x( 4 ,t ) | 9r(i »t | 4)&Zq( t ) ]  (3-115) 

for  some  4.  Thus  far,  there  are  strong  parallels  to  the  definition  of  the 
global  hypothesis  structure  in  Case  2.  However, 

PT(i,t)  *  Prob[  €>r(i  ,t  |  4)  |  ZG(  t )  )  .  (3-116) 

The  definition  of  Pj(i,t)  is  ad  hoc,  but  carefully  chosen,  as  will  be  seen. 

The  number  of  track  hypothesis  structures  at  t  is  denoted  by  Mx(t). 

To  define  how  the  components  of  track  hypothesis  structures  are  calcu¬ 
lated,  it  is  convenient  to  define  the  look  hypothesis  structure,  denoted  by 
HL(j,t).  It  has  one  component: 
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HlU.O  *  (^l(  t )  j  ,t) )  for  t>0  (3-117) 

where 

dL(t|j,t)  »  [bL(t|j,t),iL(t|j,t)]T  for  t>0  (3-118) 

and  where  dj/t  |j  ,t)  is  the  value  of  d'(£»t)  assumed  by  the  j-th  look  hypoth¬ 
esis  at  t  for  some  unspecified  l;  bj^tjj.t)  and  i  j^(  t  j  j  ,  t )  are  similarly  rela¬ 
ted  to  b(£,t)  and  i(£,t)  for  the  same  value  of  i.  The  j-th  look  hypothesis  at 

t  is  correct  if  and  only  if 

dL(t|i»t)  *  d'(tft)  for  t>0  and  some  £  .  (3-119) 

For  convenience,  we  denote  this  event  by  ©L(j,t|fc).  We  denote  the  number  of 
look  hypotheses  at  t  by  tlL(t).  On  any  scan  t,  any  one  target  can  be  "unborn," 
can  be  "undetected,"  can  be  "dead,"  or  can  cause  any  one  of  the  L(t)  detec¬ 
tions  on  that  scan.  Thus, 

ML(t)  =  L(t)  +  3  for  t>0  .  (3-120) 

Each  target  hypothesis  structure  is  constructed  from  a  sequence  of  con¬ 
secutive  look  hypothesis  structures.  Thus  each  target  hypothesis  structure  at 
t  can  be  formed  from  a  target  hypothesis  structure  at  t-1  and  a  look  hypoth¬ 
esis  structure  at  t.  We  symbolize  this  relationship  in  the  usual  manner: 

HT(i,t-l)*HL(j,t)  ->  HT(k,t)  .  (3-121) 
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As  in  Case  2,  this  generative  formula  applies  only  for  t>0  and  for  certain  i 
and  j  pairs.  For  t«0 ,  we  have  but  one  track  hypothesis  structure  Hx(l,0), 
such  that 

bT( 1,0)  =  B  (3-122) 

and 

PT(1,0)  =  1  .  (3-123) 

All  of  the  other  components  of  this  hypothesis  structure  are  undefined.  The 
indices  i  and  j  are  restricted  to  be  such  that  dj/j  ,t)  is  reachable  from 
dx(i,t-l)  in  one  step.  These  restrictions  exactly  parallel  those  in  Case  2: 

bL(t|j>t)  =  B  and  bx( t-1 | i , t-1 )  =  B 

hLCtlj.t)  *  D  and  bp( t-1 | i ,t-l )  *  X 

t>L^ t  |  j  » t >  =  D  and  b-p (t  —  l)i,t— 1)  =  D  or  D  , 


or 

bL(t|j,t)  =  X  and  bT( t-1 | i , t-1 )  -  D  or  X  .  (3-124) 

The  generative  procedure  also  closely  parallels  that  in  Case  2.  The  new 
hypothesized  event  state  sequence  ^(k.t)  is  generated  from  Dx(i,t-1)  and 
dL( j .t)  by  letting 


dp(T|k,t)  -  dx(T|i,t-l)  for  0<T<t-l 


(3-123a) 
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and 

dT(t|k,t)  -  dL(t|j,t)  .  (3-125b) 

This  definition  gives  rise  to  the  observation  that 

eT(i,t-l|£)&0L(j,t|i)  =  ©rCk.tU)  (3-126) 

for  t>0  and  £>0.  This  observation  can  be  used  to  justify  the  following  pro¬ 
cedures  for  updating  the  second  and  third  components  of  the  track  hypothesis 
structures : 

xT(t|k,t)  -  xT(t|i,t-l)+K(t|i,t-l)[z(iL(t|j,t),t)-CxT(t|i,t-l)] 

( 3- 127a) 

and 

rT(t|k,t)  =  [I-K(t|i,t-l)CT]ET(t|i,t-l)  (3— 127b) 

if  bL(t|j,t)=D  and  bj( t-1 | i ,t-l )=D  or  D, 

xT(t]k,t)  =  xT(t  ji,t-l) 

and 

Et(c |k,t)  =  EtC  1 1 i , t-1 )  , 

if  t»L.(  1 1  j  ,t)=D,  and 

xT(t|k,t)  -  K0z(lL(t|j,t),t)  (3-129a) 


( 3-128a) 


( 3— 128b) 
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IT(t|k,t)  -  [I-KoCtJZ0 


( 3— 129b) 


if  bL(t|j,t)=D  and  bj( t-1 |i ,t-l)»B  where 


xx(t|i,t-l)  =  AxX(t-l |i ,t-l)  , 


ZT( t | i , t-1 )  =  AZx(t-l |i,t-l)AT  +  W  , 


(3-130a) 


(3— 130b) 


K(t  |  i  ,  t-1 )  =  Zx(t  |i,t-l)ClV+CZT(t  Il.t-DCT]-1  ,  (3-130c) 


K0  =  IC0[V+CI0CT]-1  . 


(3-130d) 


The  update  of  PT(i,t)  depends  on  bL ( t | j , t )  and  b-p(  t-1  |  i  ,  t-1 )  as  follows: 


PT(k,t)  *  Px(i ,t-l ) 


bi_(  t  |  j  ,  t ) *B  and  bX( t-1  |  i  ,t-l )  =  B 


PX( i * t-1 ) XB[ G( z( iL(t | j ,t) ,t) ,V+CZ0CT)/pFA] 


if  bL(t|j,t)=D  and  bx( t-1 | i ,t-l )*B 


PX( i i t-1 )PD[ G(z( iL(t |j ,t),t)-CxT(t |i,t-l),V+CIT(t |i,t-l)CT)/pFA] 


if  b^(t | j ,t)“D  and  bx(t-l | i ,t-l)-D  or  D 


Px(i  ,t-l)  [  l-Pp-Px)  if  bX(  1 1  j  ,  t  )*D  and  bX  (t-1  |i,t  —  1)«=D 


PT(i,t-l)Px 


if  bj,(t  |  j  ,t)“X  and  bX( t-1 | i , t-1 )*D 


k 


PX(i,t-l)[l-PD]  if  br(t | j ,t)-D  and  bX( t-1 j i ,t-l )=D  , 
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PT(i,t-l)  if  bL(t|j,t)=X  and  bT( t-1 | i , t-1 )=X 


(3-131) 


where 


-K/2,  -1/2  -yTY_1y/2 

G(y,Y)  =  (2rr)  Y  e 


(3-132) 


The  expression  Q(t)  does  not  appear  in  this  update  unlike  previous  cases. 
Growth  in  the  number  of  target  hypothesis  structures  can  be  bounded  as 
follows : 


MT(t)  <  Mt( t- 1 ) [ L( t )+3 ]  for  t>0  . 


(3-133) 


As  with  the  global  hypothesis  structures  in  Case  2,  the  rate  of  growth  can  be 
geometric . 

The  second  step  in  this  alternative  formulation  of  the  optimal  tracking 
algorithm  requires  an  alternative  definition  of  the  global  hypothesis  struc¬ 
ture.  The  i-th  alternative  global  hypothesis  structure  has  three  components: 


H'(i,t)  -  (NG(i,t),KG(i,t),PG(i.t)) 
G 


(3-134) 


where  NG(i,t)  and  PG(i,t)  have  their  usual  meanings,  where 


Kg( i > t )  -  {<G( A | i ,t ) } 


NG(i,t) 


(3-135) 
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and  where  each  *c( £ | i , t )  designates  a  target  hypothesis  structure.  In  effect, 
each  alternative  global  hypothesis  assumes  a  number  of  targets  and  associates 
a  track,  hypothesis  with  each  one. 

There  are  two  restrictions  on  the  target  hypothesis  structure  index  set 
K<](i,t).  One  restriction  is  that 

1t( t I KG( £ | i  .t ) ,  t )  *  tf( T I <g( 1 1 »t )  ,t )  (3-136) 

for  all  i*l'  and  0<t<t.  This  restriction  simply  means  that  an  alternative 
global  hypothesis  cannot  include  two  target  hypotheses  based  on  a  common  de¬ 
tection.  The  second  restriction  is  that  for  i  and  i'  such  that  0<£<£'  <N'G(i  ,t) 
then  either 


|i»t),t)  =  B  and  bx(  t  |  kG(  £  |  i  ,  t ) ,  t )  =  D  (3-137a) 

i.e.,  one  target  is  born  before  the  otner;  or,  if  both  targets  are  born  at  the 
same  time  t 


Uz(  iT(  t  |  kG(£'  I1 1 1)  ,  t)  II  >  Jz(iT(T|iCG(^li.t)»t)11  ( 3—1 37b) 

This  restriction  simply  means  that  targets  are  indexed  first  in  the  order  of 
the  "birth"  and  second  by  the  order  of  their  "birth"  measurement's  magnitude. 
The  method  of  calculating  Pj(i,t)  was  chosen  so  that 

NG(i,t)  -  t 

PG(i,t)  -  n  PTUGU|i,t),t)  n  Q(t)  (3-138) 

4-1  t—  1 
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where 


-XB  -Afa  L(t) 

Q(t)  =  e  [e  Ipfa  /L( t) ! ] ] / Prob[ Zs( t) | Zq{ t-1 ) ] 

for  t>0  .  (3-139) 

It  remains  to  demonstrate  the  correctness  of  this  formula.  The  first  step  is 
to  note  that  the  number  of  alternative  global  hypotheses  at  t  equals  the  num¬ 
ber  of  global  hypotheses,  since  both  sets  of  hypotheses  must  include  all  pos¬ 
sible  hypothetical  associations  of  detections  with  targets.  For  convenience, 
we  assume  that  the  indices  of  the  original  and  alternative  hypothesis  struc¬ 
tures  are  equivalent.  That  is,  we  assume 

d '  ( £  , t  |  i  , t )  *  dT(  t  |  <G(  £  [  i  , t: ) ,  t  >  for  T<t  (3-140a) 

G 


and 


x'U,t|i,t)  =  1 1  kq(  £  |  i  ,t )  ,  t )  , 


( 3-1 40b) 


both  for  t>0. 

The  second  step  is  to  observe  that  Eq.  3-85b  implies 


d'U,t  j k , t )  =  d'  (fc,t  |  j  ,t) 
G  s 


(3-141) 


for  l<H<NQ(k,t)  and  for  some  j  corresponding  to  every  k.  Thus, 
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d  (£,t|j,t)  -  dT(t | <G( i |k,t ) ,t)  (3-142) 

s 

for  l<£<NG(k,t)  and  for  some  j  correspopding  to  every  k.  Using  the  relation 
between  j  and  k  implied  in  this  last  equation,  one  can  define  the  target  index 
sets  4t3b( k , t ) ,  etc.  as  in  the  last  subsection  with  a  slightly  different  inter¬ 
pretation,  e.g., 

bT(t I lk.t ) ,t)  «  D  and  br( t-1 | xG( i | k , t ) , t )  =  B  (3-143) 

if  and  only  if  £e'fDB(k,t)  • 

The  third  step  is  then  to  rewrite  Eq.  3-138  as  follows: 

PG(k.t)  =  -£  PT(KG(i|k,t),t)  ~  11  PT(Ko(i|k,t),t)  x 

ieH-DBU.t)  £efDD(k,t) 


-_n  PT(KG(i|k,t),t)  PT(KG(£|k,t),t)  x 

^ETDD(k,t)  ieTDD(k,t) 


*  £  PT(KG(i|k,t),t)  PT(KG(£|k,t),t  x 

ie'fDD(k,t)  2eH'DD(k,t) 


'  "  PT(KG(£|k,t) ,t)  n  Q(t)  (3-144) 

T=1 


Using  Eqs .  3  125  and  3—131,  it  is  possible  to  further  the  rewriting  process: 
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NG(i,t-l)  -  t-1 

PG(k ,t )  =  n  PT(.cG(m|i,t-l),t-l)  n  Q(t)  X 

m= 1  x=  1 


Ab[  G(z(  ij(  t  j  icG(  £  j  k  ,  t )  ,  t  )  ,  V+CI0CT)/ppa]  x 

*e'PDb(k,t ) 


11  _  PD[G(z(  iT(t  |  <G(  £|k,t )  ,t)-CxT(t  |  kG(£  |  i  ,t-l)t-l)  , 

,  t  )U>[)D(k  ,  t ) 


V+CIT(t|<G(£ji)t-l))C-l)CTj  x 


^DD^k.t)  to(k,t)  4X)D(k,t) 

(1-PD-PX]  PX  [1-PdI  Q(t)  (3-145) 


for  some  i  corresponding  to  each  k.  Because  the  first  two  factors  of  this 
last  equation  simply  equal  PG(i,t-l),  the  correctness  of  Eq.  3-138  is  demon¬ 
strated  intuitively. 

The  motivation  for  this  alternative  form  for  the  hypothesis  structures 
and  for  the  optimal  algorithm  is  not  obvious  at  this  point.  It  is  in  simpli¬ 
fying  and  then  implementing  the  algorithm  that  the  advantages  reveal  them¬ 
selves.  The  advantages  are  partially  illuminated  in  the  next  section  which 
sketches  an  implementation  of  the  optimal  algorithm. 
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SECTION  4 

ALGORITHM  DEVELOPMENT:  PRACTICAL  ALGORITHMS 

Defining  an  algorithm  which  is  optimal  given  all  of  the  phenomena  which 
must  be  captured  (false  and  missing  detections,  target  maneuvers,  data  associ¬ 
ation,  etc.)  is  difficult  even  in  the  absence  of  computational  constraints. 

To  take  such  an  algorithm  and  make  it  computationally  practical  without  sacri¬ 
ficing  much  performance  is  also  difficult.  Practicality  requires  that  the 
number  of  hypotheses  actually  generated  and  tested  by  a  tracking  algorithm  be 
bounded  by  some  value  dependent  on  the  host  computer.  Thus  some  hypotheses 
must  be  ignored.  This  fact  leads  to  the  necessity  of  defining  rules  for  the 
pruning,  combining,  and  ignoring  hypotheses  in  an  efficient  fashion  while 
maintaining  a  high  level  of  performance.  Here  we  must  be  careful  in  defining, 
even  informally,  what  we  mean  by  a  "high  level  of  performance."  One  obvious 
measure  of  performance  is  the  probability  that  the  correct  hypothesis  is  not 
ignored  or  discarded.  However,  if  one  ignores  or  discards  even  a  single 
hypothesis  there  is  a  nonzero  probability  that  it  will  be  the  true  hypothesis 
that  is  ignored  or  discarded.  Furthermore,  if  all  but  a  bounded  number  of 
hypotheses  are  ignored,  then  the  probability  that  the  true  hypothesis  is 
ignored  or  discarded  grows  as  t+<®,  and,  although  we  have  not  yet  attempted  to 
prove  it,  we  conjecture  that  in  many  cases  this  probability  approaches  unity 
as  t-*-®. 
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It  is  very  important  that  the  last  sentence  be  understood  thoroughly,  as 
it  clearly  points  out  the  inadequacy  of  using  the  probability  of  keeping  the 
true  hypothesis  as  the  critical  measure  of  performance  and  in  fact  suggests 
both  alternative  measures  of  performance  and  potential  pruning  strategies. 
Specifically,  in  our  Case  3,  the  true  hypothesis  includes  the  precise  specifi¬ 
cation  of  exactly  when  each  target  is  born  and  dies,  the  complete  list  of 
scans  on  which  each  target  is  detected,  and  the  data  associations  for  detected 
targets.  If  one  were  to  include  maneuvers  and  formations,  the  true  hypothesis 
would  include  the  precise  times  and  descriptions  of  every  maneuver  of  every 
target  and  the  precise  times  of  the  splitting  of  targets  from  formation  or  the 
merging  of  targets  into  formation.  With  this  definition,  an  hypothesis  cor¬ 
responding  to  absolute  truth,  except  for  a  single  misclassif ication  of  one 
return  as  a  false  alarm  rather  than  a  detected  target,  is  a  false  hypothesis, 
even  if  this  one  error  corresponds  to  an  event  several  hundred  scans  in  the 
past!  From  this  example  it  is  clear  that  there  are  many  hypotheses  that  are 
close  to  the  truth  and  consequently  will  produce  estimates  virtually  indis¬ 
tinguishable  from  that  of  the  true  hypothesis,  except,  perhaps,  for  small 
intervals  of  time  corresponding  to  points  at  which  the  "nearby"  hypothesis  is 
in  error.  Thus,  one  is  led  to  the  notion  of  the  resilience  of  a  suboptimal 
MMAE  algorithm.  There  are  a  variety  of  ways  in  which  one  can  quantify  this 
intuitive  idea.  For  example,  one  might  take  as  a  measure  of  resilience  the 
probability  that  at  least  one  hypothesis  "sufficiently  close”  to  the  true 
hypothesis  is  retained  or  the  expected  percentage  of  time  that  such  an  hypoth¬ 
esis  is  retained. 

While  we  have  not  pursued  the  detailed  calculation  of  such  measures  or 
the  optimal  choice  of  pruning  rules  (which  is  a  very  difficult  problem). 
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having  specified  the  optimal  algorithm  provides  us  with  an  excellent  founda¬ 
tion  both  for  devising  promising  rules  and  for  evaluating  their  performance, 
both  through  simulations  and  eventually  via  analysis.  To  date  our  work  has 
followed  the  forms  of  these  methods  of  evaluation,  but  before  we  turn  to 
describing  several  rules  and  subsequently  the  evaluation  of  their  performance, 
let  us  make  several  comments  about  analytical  methods.  Specifically,  to 
develop  such  methods,  what  is  needed  are  two  types  of  tools: 

1 .  methods  for  measuring  the  distance  between  hypotheses  and 
for  determining  the  set  of  hypotheses  that  are  "suffi¬ 
ciently  close  to  the  truth;"  and 

2.  methods  for  evaluating  the  probability  that  all  hypotheses 
that  are  "sufficiently  close"  are  discarded  when  particu¬ 
lar  pruning  rules  are  used. 

While  neither  of  these  tasks  is  an  easy  one,  there  are  existing  tools  that 
should  prove  of  value  for  each.  Specifically,  as  we  have  indicated,  "suffi¬ 
ciently  close"  can  be  defined  in  terms  of  estimation  accuracy.  Thus,  standard 
results  on  suboptimal  linear  filtering  can  be  used.  Alternatively,  one  might 
consider  using  the  recent  results  of  Baram  [8],[9],[18J  and  others  on  distance 
measures  between  different  models  (i.e.,  hypotheses)  and  on  the  selection  of 
suboptimal  modes  for  filter  design.  The  problem  of  evaluating  probabilities 
for  specific  pruning  rules  is  essentially  a  complex  level  crossing  problem,  as 
all  such  rules  can  be  viewed  as  consisting  of  a  set  of  threshold  or  comparison 
tests.  There  exists  a  large  literature  on  this  topic,  and  in  particular  the 
recent  computationally  efficient  approximate  methods  of  Gallager  and  Helstrom 
[19],  and  researchers  at  ALPHATECH  [20]  working  on  a  different  project  should 
be  of  value.  In  the  remainder  of  this  section  we  focus  our  attention  on 
several  specific  pruning  rules  which  have  clear  and  obvious  advantages  and 
interpretations . 
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One  type  of  rule  that  might  be  used  is  to  explicitly  discard  track, 
hypotheses  using  much  the  same  criteria  as  are  used  implicitly  in  currently 
fielded  airborne  surveillance  radar  systems.  Those  systems  often  employ 
tracking  “correlation  gates"  to  decide  which  detection  should  be  used  to  up¬ 
date  which  tracks.  In  terms  of  the  algorithm  described  in  subsection  3.2.4, 
such  gating  would  discard  the  k-th  track  hypothesis  if,  in  the  course  of 
generating  Pf(k,t)  from  Px(i,t-1)  as  per  Eqs.  3-131  and  3-132 

G(z(iL(t|j,t),t)-CxT(t|i,t-l),V+C2:T(t|i,t-l)CT)  >  TG  .  (4-1) 

In  essence,  this  equation  establishes  an  elliptical  track  correlation  gate 
about  the  predicted  measurement  coordinates  Cxf(c  ji ,t-l) ,  with  orientation  and 
axes  determined  by  the  covariance  of  that  prediction  and  of  the  measurements. 
For  this  correlation  gating  scheme,  one  can  explicitly  determine  the  probabil¬ 
ity  of  discarding  the  correct  hypothesis  on  any  particular  scan  with  this 
test;  it  is 

e_TG/2  .  (4-2) 

This  result  illustrates  the  advantage  of  using  the  optimal  algorithm  as  a 
starting  point  for  defining  a  tracking  algorithm  even  if  familiar  ad  hoc 
approaches  are  used  to  reduce  storage  and  computational  requirements.  Specif¬ 
ically,  having  the  optimal  algorithm  available  allows  one  to  ascertain  pre¬ 
cisely  the  nature  of  the  suboptimality  introduced  by  any  particular  decision 
rule  and  to  quantify  the  extent  of  suboptimality. 

The  approach  in  the  preceding  paragraph  discards  track  hypotheses.  We 
presume  every  global  hypothesis  which  involves  a  discarded  track  hypothesis  is 
also  discarded.  Another  technique  can  be  adopted  which  begins  by  discarding 
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global  hypotheses  and  then  goes  on  to  discard  track  hypotheses  which  are  not  a 
component  of  some  surviving  global  hypothesis.  In  this  case  one  retains  the 
Tm  most  probable  global  hypotheses  and  discards  all  the  rest. 

While  both  of  these  approaches  should  reduce  to  a  trivjal  level  the  prob¬ 
ability  of  discarding  the  correct  global  hypothesis  if  Tg  and  are  suffi¬ 
ciently  large,  they  still  allow  the  possibility  that  the  correct  global 
hypothesis  will  be  discarded.  However,  as  indicated  previously,  discarding 
the  correct  global  hypothesis  need  not  be  a  catastrophe  or  even  a  serious 
problem.  For  simplicity,  consider  two  track  hypotheses,  HX(l,t)  and  Hx(2,t). 
Both  are  based  on  detections  of  the  same  target,  with  one  exception.  Suppose 
that 

dx(r|l,t)  ■  dx(  1 1 2  ,  t )  ■=  d'(i,r)  (4-3) 

for  some  l  and  for  all  T<t  except  tq.  For  t  ,  let 

bX(T0|l,t)  -  bU,r0)  =  D  ,  (4-4) 

ix(T0|l,t)  =  iU,tq)  ,  (4-5) 


and 

bx( Tq | 2 ,t )  -  D  .  (4-6) 

That  is,  Hx(l,t)  is  the  correct  hypothesis  while  Hx(2,t)  differs  in  assuming 
that  the  target  was  not  detected  at  tq  although  it  actually  was. 

If  rQ“t,  then  the  state  estimates  and  associated  covariance  will  show  a 
difference,  although  not  necessarily  a  large  one  unless  the  target  was  very 
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recently  born.  As  t-Tg  becomes  larger  and  the  difference  in  the  hypotheses 
recedes  into  the  past,  the  state  estimates  and  associated  covariances  become 
closer  and  closer  in  value.  When  t-t  is  large,  the  estimates  and  covariances 
become  indistinguishable.  Thus,  if  H^(l,t)  had  been  discarded  for  some 
reason,  Hf(2,t)  would  eventually  become  essentially  equivalent.  The  point  of 
this  example  is  that  our  tracking  algorithm  is  quite  resilient;  it  can  recover 
from  the  incorrect  discarding  of  the  correct  hypothesis  with  a  transient  in 
the  state  estimates  and  associate  covariances. 

This  observation  suggests  that  storage  and  computational  resources  can  be 
saved  by  combining  hypotheses,  either  track  hypotheses  or  global  hypotheses, 
which  are  essentially  equivalent.  We  term  the  process  "hypothesis  fusion." 

The  most  obvious  criterion  for  deciding  when  to  fuse  track  hypotheses  is  a 
simple  threshold  test,  a  weighted  difference  in  state  estimates.  For  Hf(l,t) 
and  Hj(2,t),  the  hypotheses  are  fused  if 

lxT(  1 1 1 ,  t  )-xT(  1 1 2 , t )  II 2  <  TF  .  (4-7) 

[ZT(t|l,t)/2+ZT(t|2,t)/2]-1 

Then,  one  hypothesis  or  the  -her  is  discarded  and  the  other  appropriately 
modified.  Continuing  with  our  example,  a  reasonable  modification  is  to  use  as 
an  estimate  the  weighted  average  of  the  estimates  of  the  hypotheses  being 
fused  and  to  use  the  average  of  their  covariances: 

xT( 1 1 1 ,t )  <»  [ ZT( 1 1 1 ,t)/2+ZT(t 1 2  ,  t ) / 2  J  x 

[ rT( t 1 1 , t )-1 /2xT( 1 1 1 ,t)+ZT(t !2,t)-1/2xT(t |2,t)]  (4-8) 

IT(t|l,t)  <-  ZT(t  1 1,0/2  +  EtC  t  I  2  ,t  )/2  .  (4-9) 
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Also,  the  probability  mass  distributed  between  the  hypotheses  being  fused  is 
summed : 


PT(t|l,t)  <=  PT(t]l,t)  +  PTC t | 2 ,t)  ,  (4-10) 

creating  a  new  H-pO.t);  Hx(2,t)  would  be  discarded.  We  would  then  similarly 
merge  any  global  hypotheses  differing  only  by  the  two  track  hypotheses, 

HT(l,t)  and  Hx(2,t). 

An  even  simpler  "hypothesis  fusion"  technique  would  be  to  combine  Hj(l,t) 
and  Hf(2,t)  if 

df(T|l,t)  ■  dx(T|2,t)  for  t-Tx<T<t  .  (4-11) 

Then  the  two  hypotheses  would  be  merged  as  in  Eqs.  4-8  through  4-10,  again 
producing  a  new  HT(l,t)  and  discarding  the  old  Hx(2,t).  Some  global  hypoth¬ 
eses  would  then  be  merged  as  needed. 

In  judging  the  efficacy  of  any  particular  rule  for  reducing  the  number  of 
hypotheses  that  are  kept  one  must  consider  several  issues.  The  first  of  these 
is,  of  course,  the  effectiveness  of  the  rule  in  reducing  the  number  of  hypoth¬ 
eses  in  a  resilient  manner.  As  indicated  earlier  in  this  section  we  have  not 
developed  detailed  analytical  tools  for  evaluating  resilience  (although  we 
have  briefly  outlined  how  this  might  be  done).  Rather,  in  Section  6  we  will 
demonstrate  the  resilience  and  qualitative  behavior  of  algorithms  using  some 
of  these  rules  by  describing  a  series  of  simulations  that  have  been  performed. 

The  second  issue  is  the  complexity  of  the  additional  calculations  re¬ 
quired  for  the  implementation  of  the  rule.  While  we  have  not  examined  this 
issue  in  great  detail  for  any  of  the  rules  we  have  proposed,  it  is  not  diffi¬ 
cult  to  see  the  key  differences  among  these  rules  in  terms  of  complexity. 


83 


ALPHATECH,  INC. _ 

Gating  rules  as  in  Eq.  4-1  offer  the  greatest  computational  savings.  This  is 
not  only  because  the  calculations  in  Eq.  4-1  are  simple.  In  fact  the  critical 
point  is  that  this  rule  determines  hypotheses  which  need  not  be  generated  at 
all .  Such  a  rule  is  termed  a  "rule  for  ignoring,"  and  additional  rules  of 
this  type  should  be  considered  if  for  no  other  reason  than  the  computational 
savings  they  may  offer.  As  a  second  example,  consider  the  rule  which  keeps 
all  but  a  specified  number  of  the  most  probable  hypotheses.  This  is  a  "rule 
for  pruning"  which  requires  that  hypotheses  be  generated  before  being  dis¬ 
carded.  Such  rules  are  no  doubt  essential  but  intuitively  should  be  used  sub¬ 
sequent  to  using  rules  for  ignoring  to  reduce  the  number  of  hypotheses  that 
are  generated.  Finally,  the  rule  described  in  Eqs.  4-3  through  4-10  is  an 
example  of  a  "rule  for  merging"  which  is  the  most  complex  type  of  rule,  as  not 
only  does  one  generate  hypotheses  but  one  also  must  search  for  hypotheses  to 
be  merged  and  then  perform  operations  to  combine  the  hypotheses.  Again  such 
rules  will  certainly  be  necessary,  and  in  fact,  they  are  of  critical  impor¬ 
tance  in  assuming  that  pruning  rules  work  in  a  desirable  fashion. 

The  third  issue  is  a  bit  more  subtle,  as  it  involves  the  intuitive  sim¬ 
plicity  of  the  rule.  Specifically,  in  a  large-scale  application  in  which 
numerous  rules  are  used,  a  critical  element  in  understanding  an  overall  algo¬ 
rithm  is  a  precise  understanding  of  what  each  individual  rule  is  meant  to  do. 
This  is  especially  important  in  algorithm  development,  modification,  and 
troubleshooting:  there  are  too  many  existing  decision  tree  algorithms  in 

which  the  sheer  complexity  of  the  system  precludes  one  being  able  to  deduce 
how  the  algorithm  will  respond  under  any  specified  conditions;  by  starting 
from  the  complete,  optimal  algorithm  and  proceeding  systematically  we  have 
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begun  to  develop  a  set  of  rules  that  will  allow  us  to  avoid  this  difficulty 
and  to  determine  a  set  of  rules  which  together  reduce  overall  complexity  to  a 
satisfactory  level  while  maintaining  a  high  level  of  performance. 
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SECTION  5 
IMPLEMENTATION 


5.1  INTRODUCTION 

The  multi  target  tracker  has  two  major  tasks:  generation  of  data  associa¬ 
tions,  and  generation  of  consistent  interpretations  of  data  associations.  To 
accomplish  these  tasks,  the  tracker  maintains  and  processes  corresponding  data 
structures:  the  association  list,  and  the  global  hypothesis  list.  The  asso¬ 

ciation  list  stores  numerical  data  which  describes  postulated  trajectories. 

The  global  hypothesis  list  stores  combinations  of  trajectories  each  of  which 
is  given  a  numerical  likelihood.  At  any  time,  the  global  hypothesis  of 
highest  likelihood  summarizes  the  tracker's  best  estimate  of  the  number  of 
trajectories  as  well  as  the  current  location  and  velocity  of  the  targets 
traversing  these  trajectories.  The  association  list  is  updated  once  in  each 
scan  cycle.  The  global  hypothesis  list  is  generated  anew  in  each  cycle  by  ex¬ 
panding  the  list  which  was  present  at  the  end  of  the  previous  scan. 

In  a  supporting  role  is  the  trajectory  list  which  stores  the  coordinates 
of  simulated  returns  and  false  detections.  Based  on  parameters  entered  by  the 
user,  the  trajectory  generator  produces  lists  of  radar  returns  with  Gaussian 
random  noise  added.  The  trajectory  list  is  a  dually-linked  list.  Individual 
trajectories  are  generated  over  a  series  of  scans,  a  single  return  in  each 
scan.  The  same  data  is  presented  to  the  tracker,  however,  as  sets  of  returns, 
one  set  per  scan.  The  parameters  upon  which  trajectories  are  based  are  not 
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available  to  the  tracker  itself.  All  elements  of  the  description  of  trajec¬ 
tories,  including  the  number  of  trajectories,  must  be  deduced  by  the  genera¬ 
tion  and  ranking  of  global  hypotheses. 

The  operation  of  the  trajectory  generator  is  straightforward,  but  the 
processing  of  the  association  list  and  the  global  hypothesis  list  is  quite 
complex.  The  rest  of  Section  5  describes  the  building  of  hypotheses  and  the 
devices  employed  to  limit  their  number. 

5.2  THE  ASSOCIATION  LIST 

The  association  list  is  constructed  from  the  primitive  elements  referred 
to  as  track  nodes  or  track  hypotheses.  A  track  node  stores  the  information 
which  uniquely  describes  one  possible  trajectory:  the  estimated  position  and 
velocity  of  the  target,  the  estimated  errors  in  these  estimates,  and  the 
numerical  likelihood  of  this  target  state.  Procedure  ASSOCIATE,  which  pro¬ 
cesses  this  list,  calls  most  of  the  computational  submodules  of  the  tracker. 
These  include  the  Kalman  filter,  the  likelihood  calculator,  and  lower-level 
matrix  mathematics  routines. 

Each  simulated  return  is  represented  in  two  ways  in  the  association  list, 
first  of  all  as  a  new  target,  and  secondly  as  a  group  of  continuations  of 
existing  tracks.  A  track  node  which  represents  a  new  target  becomes  the  root 
node  of  a  potential  tree  of  track  nodes  -  a  track  tree.  As  track  continua¬ 
tions,  the  return  is  multiply  represented  as  track  nodes  which  occupy  new 
branches  on  preexisting  track  trees.  Each  terminal  track  node  in  each  track 
tree  receives  one  additional  branch  for  each  simulated  return  (an  exception  to 
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this  is  noted  below).  The  Kalman  filter  bases  its  state  estimate  for  a  new 
branch  node  on  the  updated  state  estimate  associated  with  the  parent  node  and 
the  coordinates  of  the  return. 

Besides  the  branch  nodes  associated  with  returns,  two  additional  branch 
nodes  are  added  to  the  descendants  of  some  parent  track  nodes,  one  node  which 
corresponds  to  a  postulated  missed  detection,  and  a  second  which  corresponds 
to  a  track  termination.  Parent  nodes  representing  missed  detections,  and 
nodes  corresponding  to  track  terminations  are  given  special  handling.  Missed 
detections  are  not  given  branches  corresponding  to  terminations.  Terminations 
are  given  only  terminations;  they  never  receive  branches  corresponding  to 
returns  or  to  missed  detections. 

The  number  of  track  trees  at  the  end  of  any  scan  is  equal  to  the  total 
number  of  returns  up  to  that  time.  The  depth  of  a  track  tree  equals  the  num¬ 
ber  of  scans  which  have  been  processed  beginning  with  the  scan  in  which  the 
root  node  entered  the  association  list  up  until  the  present  scan.  The  breadth 
of  a  track  tree  is  a  function  of  its  depth,  and  the  number  of  returns  per 
scan.  Figure  5-1  traces  the  evolution  of  the  association  list  for  a  simple 
case . 


5.3  THE  GLOBAL  HYPOTHESIS  LIST 

The  global  hypothesis  list  is  a  bilevel  list.  Each  entry  in  the  list,  a 
global  hypothesis,  is  itself  a  list  of  elements  each  of  which  stores  the 
address  of  a  track  node.  These  elements  are  referred  to  here  as  hypothesis 
nodes.  The  entire  global  hypothesis  list  stores  these  lower-order  lists  in 
descending  order  of  hypothesis  likelihood. 
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Figure  5-1.  Evolution  of  the  Association  List. 
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The  Crack  nodes  to  which  hypothesis  nodes  refer  are  terminal  track  nodes 
on  some  branch  of  a  track  tree.  In  any  global  hypothesis,  each  hypothesis 
node  selects  a  track  node  in  a  unique  track  tree.  Since  each  track  tree  cor¬ 
responds  to  a  different  postulated  target,  each  hypothesis  node  refers  to  a 
different  postulated  trajectory.  If  a  global  hypothesis  postulated  N  trajec¬ 
tories,  it  will  contain  N  hypothesis  nodes  referring  to  track  nodes  in  N  dif¬ 
ferent  track  trees.  The  likelihood  of  a  track  hypothesis  equals  the  sum  of 
the  likelihoods  assigned  to  each  track  node  in  the  track.  The  likelihood  of  a 
global  hypothesis  equals  the  sum  of  the  likelihoods  of  the  track  hypotheses  to 
which  it  refers.  (Throughout  this  discussion,  the  term  "likelihood”  should  be 
interpreted  as  the  natural  logarithm  of  the  likelihood  value.) 

Figure  5-2  illustrates  the  appearance  of  the  global  hypothesis  list  after 
two  scans  with  a  single  return  in  each  scan.  No  attempt  is  made  to  represent 
ranking  by  likelihood.  Note  that  global  hypotheses  G2  and  G3  in  Fig.  5-2  are 
present  in  two  versions,  one  which  does  not  postulate  the  presence  of  a  new 
target  (G2,G3)  and  a  second  version  which  does  (G5,G6).  Each  global  hypoth¬ 
esis  is  produced  as  multiple  version  equal  in  number  to  the  number  of  combina¬ 
tions  of  possible  new  targets,  whenever  these  versions  are  not  inconsistent 
with  the  total  number  of  returns  in  the  scan.  Referring  again  to  Fig.  5-2, 
you  will  note  that  there  is  no  global  hypothesis  which  include  references  to 
both  the  D  node  of  track  tree  1  and  the  B  node  of  track  tree  2.  This  hypoth¬ 
esis,  which  states  that  the  return  is  both  a  continuation  of  a  previous  track 
and  a  new  target,  is  logically  inconsistent.  Constraints  built  into  the 
global  hypothesis  generator  (logical  constraints  -  discussed  in  subsection 
5.4)  screen  out  logically  inconsistent  hypotheses. 
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Figure  5-2.  Global  Hypothesis  List  After  Two  Scans 
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The  trajectories  postulated  by  the  six  global  hypotheses  in  Fig.  5-2  and 
their  interpretations  are  given  as  follows: 

Gl.BD:  one  target,  born  in  scan  1,  observed  in  scan  2; 

G2.BN:  one  target,  born  in  scan  1,  not  observed  in  scan  2; 

G3.BX:  one  target,  born  in  scan  1,  no  longer  existing  at 

scan  2; 

G4.B:  one  target,  seen  first  in  scan  2; 

G5.BNB:  two  targets;  one  born  in  scan  1  and  missed  in 

scan  2;  one  born  in  scan  2;  and 

G6.BXB:  two  targets;  one  born  in  scan  1,  now  dead;  one  born 

in  scan  2. 

These  hypotheses  implicitly  postulate  different  numbers  of  false  alarms: 

Gl.  0  false  alarms; 

G2.  1  false  alarm,  in  scan  2; 

G3.  1  false  alarm,  in  scan  2; 

G4.  1  false  alarm,  in  scan  J; 

G5.  0  false  alarms;  and 

G6.  0  false  alarms. 

Figure  5-3  illustrates  how  one  of  these  hypotheses,  Gl  in  Fig.  5-2,  is 
expanded  into  a  family  of  global  hypotheses  in  the  next  scan.  Note  that  no 
descendant  of  Gl  makes  reference  to  track  tree  2;  none  of  its  descendants 
postulates  any  trajectory  which  begins  in  scan  2.  During  expansion  of  G5  or 
G6,  however,  this  would  not  be  the  case. 

Gll.BDD:  one  target,  observed  in  three  scans; 

G12.BDN:  one  target,  observed  in  scans  1  and  2; 
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NOTE:  G1  IS  THE  PARENT  GLOBAL  HYPOTHESIS. 
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Figure  5-3.  Expansion  of  a  Global  Hypothesis 
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G13.BDX:  one  target,  observed  for  two  scans,  but  now  dead; 

G14.BDXB:  one  target  appearing  anew  in  scan  3;  a  second  target 
observed  in  scans  1  and  2,  now  dead;  and 

G15.BDNB:  two  targets,  one  observed  in  scans  1  and  2,  a  second 
appearing  first  in  scan  3. 

The  generation  of  global  hypotheses  is  the  responsibility  of  procedure 
PARTITION  which  calls  the  recursive  procedure  SAMPLE.  SAMPLE  behaves  as  if  it 
is  constructed  from  nested  loops,  one  loop  for  each  track,  tree  to  which  a 
parent  hypothesis  refers.  The  loops  pass  through  all  track  nodes  which  are 
descendants  of  the  track  node  in  each  tree,  which  belongs  to  the  set  of  track 
nodes  referenced  by  the  parent  hypothesis.  The  inner-most  loop  accesses  the 
highest  numbered  track  tree  and  varies  most  rapidly.  In  creating  a  descendant 
hypothesis,  one  track  node  is  selected  from  each  tree.  The  process  is 
repeated  until  all  relevant  track  nodes  in  all  legal  combinations  are  incor¬ 
porated  into  descendant  hypotheses.  This  technique  constitutes  a  depth-first 
search  of  the  association  list. 

5.4  IMPLEMENTATION  OF  CONSTRAINTS 

The  "optimal"  tracker  considers  interpretations  of  radar  returns  exhaus¬ 
tively  in  the  fashion  outlined  in  subsections  5.2  and  5.3.  Unless  heuristic 
constraints  are  selected  by  the  user,  every  internally  consistent  data  inter¬ 
pretation  is  generated,  and  assigned  a  numerical  likelihood  value.  For  all 
but  the  most  trivial  cases,  the  operation  of  the  tracker  in  optimal  mode  has 
proven  to  be  a  near  impossibility  owing  to  the  demands  on  memory  and  CPU  time. 
When  heuristics  are  applied,  the  tracker  runs  in  real  time  in  many  cases. 

Table  5-1  documents  the  rapid  growth  in  the  number  of  hypotheses  with  time, 
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and  the  effectiveness  of  applying  only  a  portion  of  the  constraints  available 
to  the  tracker,  in  limiting  the  growth  of  global  hypotheses. 

TABLE  5-1.  GROWTH  OF  THE  GLOBAL  HYPOTHESIS  LIST 
CASE:  ONE  SIMULATED  TRAJECTORY;  ONE  FALSE  ALARM  PER  SCAN 


Scan  Number 


Global  Hypotheses 


Without  Constraints  With  Constraints* 


1 

4 

(2) 

4 

(2) 

2 

62 

(10) 

6 

(3) 

3 

1,620 

(34) 

16 

(6) 

4 

- 

40 

(12) 

5 

- 

48 

(22) 

*The  constraints  include  gating,  elimination  of  track 
initializations  not  followed  by  a  detection,  and  elimi¬ 
nation  of  postulated  trajectories  which  contain  more 
than  two  successive  missed  detections.  The  number  of 
track  hypotheses  is  shown  in  parentheses. 


Even  the  optimal  tracker,  which  operates  without  heuristic  constraints, 
must  apply  logical  constraints  on  the  hypotheses  building  in  order  to  avoid 
logically  inconsistent  hypotheses.  There  are  three  logical  constraints  which 
are  applied  to  hypotheses  building,  they  are  as  follows. 

1.  No  hypothesis  may  postulate  a  greater  number  of  targets 
than  the  total  number  of  returns  in  that  scan.  Returns 
are  counted  and  a  running  total  of  the  observed  postulated 
by  a  global  hypothesis  under  construction  is  compared 
against  this  count. 

2.  No  hypothesis  may  postulate  a  greater  number  of  trajec¬ 
tories  than  the  total  number  of  past  and  present  returns. 

The  fact  that  the  number  of  track  trees  cannot  exceed  the 
total  number  of  returns  ensures  that  this  constraint  is 
obeyed. 
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3.  No  returns  may  be  assigned  to  more  than  one  trajectory 

postulated  by  one  global  hypothesis.  To  prevent  this  type 
of  inconsistency,  the  tracker  assigns  arbitrary  indices  to 
the  returns  in  a  scan.  Once  a  return  is  assigned  to  a 
postulated  trajectory,  the  index  of  the  return  is  added  to 
a  set  of  disallowed  indices.  The  set  is  examined  each 
time  a  return  is  considered  when  building  a  global 
hypothesis. 

Logical  constraints  are  not  sufficient  to  prevent  the  unmanageable  growth 
in  the  number  of  hypotheses.  Additional  constraints  -  heuristic  constraints  - 
must  be  added  to  limit  the  size  and  complexity  of  the  association  list,  and 
the  number  of  global  hypotheses.  The  choice  of  heuristic  constraints  is  arbi¬ 
trary,  but  the  following  seem  reasonable  and  have  proven  effective. 

1.  Gating:  Prior  to  adding  to  the  association  list  a  branch 
which  represents  a  detected  target,  a  numerical  parameter 
is  calculated  which  reflects  the  plausibility  of  the  esti¬ 
mated  target  position  relative  to  its  estimated  position 
in  the  previous  scan.  This  gating  parameter,  which  is  a 
function  of  the  filter  innovation  and  estimated  error  in 
this  innovation,  is  compared  to  a  preselected  constant. 

If  the  gating  parameter  exceeds  the  gating  constant,  the 
data  association  is  rejected.  The  magnitude  of  the  gating 
constant  controls  the  rigor  of  the  gating  constraint. 

Gating  is  an  especially  powerful  constraint  because  it  in¬ 
tervenes  at  a  stage  prior  to  global  hypothesis  building. 

Elimination  of  a  single  track  node  results  in  the  elimina¬ 
tion  of  numerous  global  hypotheses.  For  tracking  multiple 
targets,  gating  is  most  effective  with  trajectories  which 
are  widely  separated  relative  to  the  error  in  the  measure¬ 
ments.  In  this  case,  only  a  minority  of  potential  data 
associations  pass  the  gate. 

2.  Born  Screening:  Recall  that  the  tracker  initiates  a  new 
track  tree  for  each  return  In  each  scan.  Under  this  con¬ 
straint,  the  tracker  will  reject  the  root  node  of  this 
track  tree  on  the  following  scan  unless  there  is  a  return 
which  forms  a  gate-passing  data  association  for  the  postu¬ 
lated  target.  The  tracker  requires  confirmatory  evidence 
in  order  to  maintain  a  track  initialization.  This  con¬ 
straint  has  a  strong  effect  on  the  complexity  of  the  asso¬ 
ciation  list  by  eliminating  tracks  at  an  early  stage  in 
their  existence. 
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3.  Limitation  on  Missed  Detections:  Under  this  constraint, 
any  postulated  trajectory  can  be  given  no  more  than  a  pre¬ 
selected  number  of  missed  detections  as  continuations. 

This  constraint,  which  has  the  effect  of  removing 
spuriously  postulated  trajectories,  is  most  effective  when 
the  false  alarm  rate  is  high.  Under  this  condition  there 
is  an  elevated  likelihood  that  incorrect  data  associations 
will  lead  to  faulty  track  initializations  followed  by 
gate-passing  continuations  which  are  also  incorrect. 
Eventually  these  tracks  will  be  eliminated  when  there  are 
no  data  associations  to  support  them. 

4.  Pruning:  This  constraint  is  a  direct  limitation  on  the 
number  of  global  hypotheses.  Only  a  preselected  number  of 
hypotheses  are  retained  at  the  end  of  each  scan:  those 
hypotheses  which  are  highest  in  likelihood.  Pruning  has  a 
retroactive  effect  on  the  complexity  of  the  association 
list;  any  terminal  track  node  which  is  referred  to  by  no 
surviving  global  hypothesis  is  ignored  when  the  associa¬ 
tion  list  is  expanded  in  the  next  scan.  To  implement 
retroactive  pruning,  the  tracker  keeps  a  count  in  each 
track  node  cf  the  number  of  global  hypotheses  which  refer¬ 
ence  it.  When  a  global  hypothesis  is  discarded  during 
pruning ,  the  count  in  each  of  its  referenced  track  nodes 
is  decremented  by  one.  Only  track  nodes  with  counts  ex¬ 
ceeding  zero  are  expanded  further. 


5.5  SUMMARY  AND  CONCLUSIONS 

The  multitarget  tracker  must  form  data  associations  and  multiple  inter¬ 
pretations  of  data  associations,  employing  what  might  be  termed  a  deep  analy¬ 
sis.  This  means  that  the  interpretation  of  events  in  a  given  scan  is  heavily 
determined  by  interpretations  developed  in  previous  scans.  One  of  the  prin¬ 
cipal  tasks  of  the  program,  therefore,  is  to  maintain  a  data  structure  which 
is  capable  of  representing  the  relationships  among  observations  which  are  dis 
tributed  over  scan  time.  The  association  list  provides  the  necessary  struc¬ 
turing  of  raw  data  and  results  of  computations. 

In  developing  interpretations  of  the  information  stored  in  the  associa¬ 


tion  list,  that  is  to  say,  the  building  of  global  hypotheses,  the  program 
carries  out  what  is  essentially  a  depth-first  search.  The  efficiency  of  this 
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search  is  wholly  determined  by  the  depth  and  breadth  of  the  association  list. 
The  complexity  of  this  data  structure  is  determined,  in  turn,  by  the  effec¬ 
tiveness  of  heuristic  and  logical  constraints  in  limiting  its  rate  of  growth. 
Whereas  brute  force  pruning  of  the  global  hypothesis  list  would  seem  to  be  a 
most  powerful  constraint,  in  reality  it  is  less  effective  than  heuristics 
which  focus  on  the  association  list  itself.  Pruning  determines  which  global 
hypotheses  are  saved .  Gating  and  associated  constraints,  which  intervene 
before  hypothesis  building,  determine  which  hypotheses  will  be  formed . 

The  constraints  available  to  the  current  version  of  the  tracker  are  suf¬ 
ficiently  powerful  to  allow  real-time  performance  except  in  especially 
challenging  cases:  densely  packed,  crossing  trajectories.  For  more  complex 
models  of  the  radar  environment  and  target  behavior,  however,  additional 
heuristic  constraints  will  be  necessary.  These  will  most  likely  involve  more 
aggressive  restriction  on  the  expansion  of  the  association  list,  and  addi¬ 
tional  postprocessing  of  the  global  hypothesis  list.  In  addition  to  pruning, 
the  merging  of  nearly  identical  global  hypotheses  should  be  valuable.  In 
practice,  the  tracker  has  tended  to  produce  and  retain  families  of  hypotheses 
which  postulate  nearly  identical  target  coordinates.  Reducing  these  families 
to  single  hypotheses  should  not  impair  performance  and  will  lead  to  a  greatly 
simplified  association  list. 
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SECTION  6 

PERFORMANCE  EVALUATION 

In  this  section  we  examine  the  performance  of  the  surveillance  and  track¬ 
ing  algorithm.  In  order  to  demonstrate  the  algorithm's  major  attributes,  we 
select  four  representative  scenarios,  based  on  an  ascending  order  of  difficul¬ 
ty.  Each  of  these  scenarios  is  simulated,  and  in  result,  the  surveillance  and 
tracking  performance  of  the  algorithm  is  established.  Common  to  all  four 
scenarios  is  the  surveillance  volume  S,  defined  as  a  20x20  kilometer  plane  at 
a  distance  of  200  kilometers  away  from  the  surveillance  radar.  For  conve¬ 
nience  we  define  the  center  of  the  surveillance  range  as  the  origin  of  a  two- 
dimensional  inertial  coordinate  system.  Thus  E  may  be  defined  as 

E  =  {x,y |-10 ,000<x<10 ,000 ,-10 ,00CKy<10 ,000}  ,  (6-1) 

where  the  positive  y-axis  direction  points  to  the  "north.” 

We  further  assume  that  the  surveillance  range  is  scanned  every  10 
seconds.  For  each  scenario  we  simulate  the  total  of  nine  scans.  The  simula¬ 
tion  program  generates  at  each  scan 

1.  noisy  position  (x,y)  returns  originating  from  existing 
target(s) ; 

2.  false  alarms;  i.e.,  an  arbitrary  number  of  background 
noise  returns,  each  of  which  is  independent  and  uniformly 
distributed  over  the  surveillance  range;  and 

3.  misses;  i.e.,  the  hypothesized  target(s)  does  not  produce 
a  return.  The  accumulated  number  of  misses  for  the  entire 
trajectory  of  a  particular  target  is  arbitrary. 
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6.1  SCENARIOS 

In  the  following  we  list  the  details  of  the  four  scenarios  examined.  In 
all  cases  it  is  assumed  that  the  target(s)  velocity  is  200  m/s. 

1.  Two  crossing  targets,  no  false  alarms,  no  misses:  the 
location  of  the  targets  at  the  first  scan  is  (-6500,  6500) 
and  (-6500,  -6500)  and  their  headings  are  135°  and  45°, 
respectively .  The  heading  angle  is  measured  from  the 
“north''  direction. 

2.  Two  crossing  targets,  two  false  alarms/scan  in  first  eight 
scans,  no  misses:  the  initial  coordinates  of  the  targets 
are  the  same  as  in  Scenario  1. 

3.  Three  targets,  one  false  alarm/scan,  some  misses:  two  of 
the  targets  run  in  parallel,  and  are  separated  by  3000 
meters.  The  third  target  crosses  the  two  parallel  trajec¬ 
tories.  Target  #1  (number  assigned  arbitrarily)  is  missed 
at  the  fifth  scan,  target  #2  at  the  third,  and  target  #3 
at  the  seventh.  The  initial  coordinates  and  headings 
(i.e.,  at  the  first  scan)  are: 

target  ill:  (-7900,  3500),  105°; 

target  it 2:  (-7900,  500),  105°;  and 

target  #3:  (-7900,  -3500),  60°. 

4.  Three  targets,  no  false  alarms,  no  misses:  Similar  to 
Case  3,  this  scenario  is  characterized  by  two  targets 
moving  in  parallel  and  a  third  one  crossing.  The  diffi¬ 
culty  in  the  present  scenario  is  the  proximity  of  the  two 
parallel-moving  targets:  we  assume  a  separation  of  1000 
meters.  The  initial  coordinates  and  headings  of  the  tar¬ 
gets  (numbered  arbitrarily)  are: 

target  it l:  (-7900,  500),  105°; 

target  It 2:  (-7900,  -500),  105°;  and 

target  #3:  (-7900,  -3500),  60°. 

Finally,  we  assume  that  the  target-originating  returns  are  normally  dis¬ 
tributed,  i.e.,  the  (x,y)  target  position  is  actually  measured  (simulated)  as 
(x+vx,y+Vy),  where 
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Vx  ~  N(0,o2)  ;  V,.  -  N(0,o2)  , 

x  J  y 


an-  we  further  assume 


(6-2) 


ox  =  5  m  ,  Oy  =  500  m  (6-3) 

This  disparity  in  the  measurement  accuracy  between  the  two  axes  emulates  a 
realistic  situation  in  which  the  accurate  measurement  is  the  range  and  the  in¬ 
accurate  is  the  line-of -sight  angle. 

6.2  ALGORITHM  PARAMETERS 

Next  we  list  the  algorithm/Kalman  filter  parameters  that  were  used  in  the 
simulation.  These  were  kept  unchanged  in  all  four  scenarios,  and  are  as 
follows : 

PFA  =  false  alarm  density  =  10-10/meter2; 

«  expected  number  of  born  targets  per  scan  =  1 0—  3 ; 

Pp  =  probability  of  detection  =  0.998;  and 
Px  «  probability  of  target  mortality  =  10-4. 

The  Kalman  filter  parameters  are  chosen  as  follows: 

Process  noise 

covariance  -  W  =  diag(0,0(10  m/sec)2, (L0  m/sec)2]; 

Measurement  noise 

covariance  «  V  =  diag[(5  m)z,(500  m)2);  and 

Initial  state 

estimate  covariance  -  Z0  =  diag[(5  m)2,(500  m)2,(200  m/s)2, (200  m/s)2]. 
The  initial  state  estimate  position  variables  are  taken  as  the  (x,y)  returns 
from  scan  #1;  the  initial  velocity  estimates  are  assumed  to  equal  zero. 
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As  indicated  in  Section  5  only  a  suboptimal  implementation  of  the  sur¬ 
veillance/tracking  algorithm  is  feasible.  We  therefore  employ  four  heuristic 
modifications  to  the  full-blown  optimal  case. 

1.  Gating:  The  gating  parameter  used  is  10,  implying  that  no 
more  than  approximately  three  standard  deviations  of  the 
innovation  process  are  allowed  when  associating  new  mea¬ 
surements  with  existing  track  hypotheses. 

2.  Track  Initialization:  The  algorithm  will  not  associate  a 
miss  or  a  track  termination  following  a  target-birth 
hypothesis . 

3.  Hisses:  If,  on  a  scan  following  two  successive  misses 
pertaining  to  a  particular  track  hypothesis,  no  return 
passes  the  associated  "gate,"  that  track  is  discounted. 

4.  Pruning:  No  more  than  10  top  ranking  (i.e.,  with  highest 
likelihood)  global  hypotheses  are  retained  after  the  pro¬ 
cessing  of  each  scan. 

6.3  PERFORMANCE  EVALUATION 

In  the  performance  analysis  we  present  figures  which  illustrate  the  algo¬ 
rithm's  surveillance  and  tracking  for  each  scenario.  In  these  figures,  the 
following  symbology  is  used. 

Returns:  O;  in  the  cases  which  involve  false 

alarms,  the  returns  are  numbered  by 
the  order  of  appearance  (scan  number); 
the  only  exception  is  Fig.  6-18  in 
which  all  returns  are  enumerated,  as 
explained  in  the  sequel; 

Hypothesized  Trajectory  a ,  4  ,  or  0  (each  symbol  pertains 

Position  Estimates:  to  a  particular  track  hypothesis); 

Estimated  Trajectories:  ....... _ ,  or  and 

Target  Trajectories:  Solid  Lines.  Targets  are  ordered  by 

encircled  numbers  (assigned  arbi¬ 
trarily)  next  to  their  trajectories. 
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6.3.1  Scenario  1 

Figures  6-1  through  6-5  Illustrate  the  surveillance/tracking  performance 
for  the  simplest  case  considered.  Shown  are  the  five  top-ranking  hypotheses, 

1. e.,  those  with  the  highest  likelihoods  of  being  correct,  after  the  ninth 
scan.  Figure  6-1  illustrates  the  Rank  1  hypothesis.  Fig.  6-2  illustrates  Rank 

2,  etc.  In  the  case  at  hand,  the  correct  hypot.  .'is  is  ranked  number  1  which 
is  the  desired  result;  i.e.,  all  returns  are  associated  by  the  algorithm,  with 
the  correct  target  positions  as  shown  in  Fig.  6-1. 

Figure  6-2  illustrates  the  Rank  2  hypothesis,  and  it  is  evident  that  at 
the  vicinity  of  the  targets'  crossing  (sixth  scan)  an  erroneous  return/target 
association  is  hypothesized.  In  other  words,  after  the  ninth  scan  the  Rank  2 
hypothesis  associates  the  sixth  scan  return  from  target  #1  with  target  #2,  and 
vice  versa. 

Figure  6-3  (hypothesis  Rank  3)  illustrates  a  different  type  of  error: 
the  Scan  1  return  which  originates  from  target  #2  is  Interpreted  as  a  false 
alarm;  consequently,  this  target  Is  assumed  to  be  born  only  on  the  second 
scan.  Naturally,  this  misinterpretation  results  in  a  different  track  estimate 
(compare  with  Fig.  6-1). 

Figure  6-4  illustrates  the  Rank  4  hypothesis.  This  hypothesis  is  similar 
to  the  Rank  2  hypothesis  but  the  erroneous  target/return  association  occurs  at 
the  fifth  rather  than  at  the  sixth  scan. 

Finally,  the  hypothesis  ranked  as  fifth  (Fig.  6-5),  assumes  that  the 
ninth  scan  return  from  target  # 2  is  a  false  alarm.  Therefore,  the  hypothesis 
is  that,  on  the  ninth  scan,  target  #2  is  not  detected. 
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Figure  6-3.  Scenario  1,  Scan  9,  Rank  3  Hypothesis. 
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Figure  6-4.  Scenario  1,  Scan  9,  Rank  4  Hypothesis. 
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In  summary,  we  have  examined  a  simple  but  an  interesting  scenario.  The 
truth  hypothesis  was  correctly  ranked  first,  and  the  subsequent  rank/hypoth¬ 
eses  illustrate  a  variety  of  plausible  track  estimates,  which,  strictly  speak¬ 
ing,  are  erroneous.  However,  note  that  the  erroneous  hypothesized  tracks  are 
relatively  similar  to  the  correct  ones. 

6.3.2  Scenario  2 

Scenario  2  is  similar  to  Scenario  1  in  the  sense  that  the  simulated  tar¬ 
gets  and  their  associated  returns  are  identical.  Presently,  however,  we  com¬ 
plicate  the  hypotheses  tree  by  augmenting  the  correct  returns  with  simulated 
false  alarms  -  two  per  scan,  in  the  first  eight  scans.  We  will  see  that 
adding  these  false  alarms  does  not  degrade  the  algorithm's  performance,  but  it 
produces  a  different  type  of  erroneous  hypotheses:  it  generates  hypothetical 
new  tracks,  which,  obviously,  are  false. 

Figure  6-6  illustrates  the  Rank  1  hypothesis  which  is  similar  to  the  Rank 
1  hypothesis  of  Scenario  1  (the  truth  hypothesis).  The  false  alarms  in  the 
first  eight  scans  (indicated  by  the  scan  number)  have  no  effect  in  this  case. 
Similarly,  the  second  ranking  hypotheses  of  the  two  scenarios  are  exactly  the 
same.  (Rank  2  hypothesis  of  the  present  scenario  is  not  shown.)  More  inter¬ 
esting  cases/hypotheses ,  however,  are  shown  in  Figs.  6-7  and  6-8. 

Figure  6-7  illustrates  the  Rank  3  hypothesis.  It  is  shown  that  targets 
#1  and  #2  are  tracked  correctly,  similar  to  the  Rank  1  hypothesis.  However,  a 
false  alarm  at  the  sixth  scan  is  taken  by  the  algorithm  as  a  birth  of  a  new 
target.  Subsequently,  this  "target"  is  "detected”  at  the  seventh  scan  by 
hypothesizing  an  association  with  another  false  return.  Note  that  these  two 
subsequent  false  returns  are  sufficiently  close  so  that  an  hypothesis 


109 


ALPHATECH,  INC 


i 


km 


10 


n4 


8  > 

6  • 

4  ■ 

2  • 


-10  L 


NINTH  SCAN  "MISS" 


o5 


o5 

o2  o4 

08  °2 


-10  -8  -6  -4  -2  0  2  4  6  8 


Figure  6-7.  Scenario  2,  Scan  9,  Rank  3  Hypothesis. 


ALPHATECH,  INC 


Figure  6-8.  Scenario  2,  Scan  9,  Rank  4  Hypothesis 


T 


ALPHATECH,  INC. _ 

associating  the  two  cannot  be  "gated-out."  On  the  eight  and  ninth  scans,  how 
ever,  this  track  makes  no  association  which  would  pass  the  gate  with  any  re¬ 
turn.  It  is,  therefore,  postulated  that  the  "target”  was  missed  on  these 
scans.  Note  that  unless  a  new  false  return  occurs  on  the  path  of  this  hypo¬ 
thetical  track  in  the  tenth  scan,  that  track  -  on  the  basis  of  our  suboptimal 
ity  assumption  -  will  be  dropped  entirely. 

Figure  6-8  illustrates  a  similar  situation  in  the  Rank  4  hypothesis.  A 
fictitious  target  is  hypothesized  to  have  been  born  by  a  false  return  on  the 
eighth  scan.  This  "target"  is  then  "detected"  in  the  ninth  scan  by  associa¬ 
ting  with  a  return  which  actually  originates  from  target  #2.  Therefore,  the 
track  associated  with  target  #2  must  assume  that,  in  the  ninth  scan,  a  missed 
detection  has  taken  place  (see  Fig.  6-5).  Note  that  all  of  the  falsely  hypo¬ 
thesized  trajectories  are  “young"  and  relatively  short  in  duration. 

6.3.3  Scenario  3 

In  this  scenario  we  assume  an  engagement  of  three  targets,  two  of  which 
move  in  parallel  and  one  of  which  is  crossing.  We  also  inject  one  false  re¬ 
turn  per  scan  for  the  first  eight  scans,  and  simulate  misses  on  the  fifth, 
third,  and  seventh  scans  for  the  first,  second,  and  third  target,  respective¬ 
ly.  The  scenario  is  not  particularly  difficult  as  the  target  returns  are 
relatively  sparse.  We  use  it,  however,  to  illustrate  several  interesting 
points : 

1.  we  illustrate  the  difficulty  in  attempting  to  perform  the 
surveillance  and  tracking  task  manually; 

2.  we  show  that  if  a  finite  number  of  scans  is  considered, 
the  truth  hypothesis  does  not  necessarily  rank  first; 
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3.  the  Rank  1  hypothesis,  albeit  erroneous,  is  not  far  from 
the  truth:  after  a  "reasonable"  number  of  scans  (greater 
than  four),  the  tracking  performance  is  adequate;  and 

4.  we  illustrate  the  mechanics  of  propagating  the  track- 
hypotheses  throughout  several  scans. 

We  start  by  inspecting  Fig.  6-9  which  illustrates  all  the  returns  from 
the  initial  eight  scans.  Clearly,  given  this  scenario,  manual  tracking  would 
be  very  difficult,  if  not  outright  impossible.  Indeed,  certain  patterns  can 
be  identified,  but  a  mental  hypothesis  ranking  seems  to  be  unattainable. 

Next,  we  examine  Fig.  6-10  which  is  the  Rank  1  hypothesis  in  the  present 
scenario.  Although  it  ranks  first,  this  hypothesis  is  not  the  truth  hypoth¬ 
esis,  as  it  fails  to  interpret  the  absence  of  a  return  from  target  It 2  on  the 
third  scan  as  a  miss.  As  a  result,  it  postulates  a  birth  of  target  It 2  only  on 
the  fourth  scan.  However,  the  tracking  performance  from  the  fifth  scan  on  is 
correct.  In  particular,  the  two  latter  misses  occurring  in  the  two  other 
trajectories  are  identified  correctly. 

Finally,  we  illustrate  the  propagation  of  the  Rank  1  hypothesis  through 
scan  numbers  2,  3,  4,  and  5  (Figs.  6-11,  6-12,  6-13,  and  6-14,  respectively). 
Only  target  It 3  is  tracked  correctly  throughout  all  four  scans.  Target  #1  is 
detected  only  on  the  third  scan,  and  target  It 2  only  on  the  fifth  (it  is 
assumed  to  have  been  born  on  the  fourth).  It  is  also  evident  that  on  the 
second  and  the  fourth  scans,  hypothetical  tracks,  effected  by  false  returns, 
are  postulated.  These,  however,  disappear  from  the  Rank  1  hypothesis 
thereafter . 

6.3.4  Scenario  4 

The  final  example  illustrates  perhaps  the  most  interesting  case,  in  which 
two  very  close  targets  move  in  parallel,  and  their  trajectories  are 
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Figure  6-9.  Eight  Scan  Returns  of  Scenario  3. 
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Scan  3,  Rank  1  Hypothesis. 
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Figure  6-13.  Scenario  3,  Scan  4,  Rank  1  Hypothesis. 
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Figure  6-14.  Scenario  3,  Scan  5,  Rank  1  Hypothesis 
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intersected  by  a  third  one.  Three  hypotheses  are  presented  in  Figs.  6-15 
through  6-17.  Figure  6-16  shows  the  truth  hypothesis  which,  in  the  case  at 
hand,  ranks  third.  The  Rank  1  hypothesis  is  shown  in  Fig.  6-15  which  is  prac¬ 
tically  identical  to  the  truth  hypothesis;  only  a  minor  crossing  between  the 
returns  from  target  #1  and  02  occurs  on  the  first  scan.  The  most  interesting 
case,  however,  is  shown  in  Fig.  6-17,  in  which  a  total  crossover  between  the 
two  parallel  targets  occurs  on  the  fourth  scan.  This  hypothesis  is  ranked 
fifth.  Clearly,  with  such  congested  returns  a  crossover  of  this  type  cannot 
be  avoided. 

6.4  PERFORMANCE  OF  THE  GATING  LOGIC 

The  final  issue  to  be  discussed  in  this  section  is  the  effectiveness  of 
the  gating  algorithm.  For  the  purpose  of  illustrating  the  performance  of  the 
gating  logic  consider  Fig.  6-18,  in  which  a  simple  tracking  case,  on  the  third 
scan,  is  shown.  We  assume  a  single  target  starting  at  (-400,  0)  and  bearing 
90°.  We  assume  no  missed  returns  but  we  simulate  a  uniformly  distributed 
false  return  on  the  second  and  the  third  scans.  All  five  returns  are  enumer¬ 
ated  as  shown.  Next,  examine  Table  6-1  which  lists  all  the  possible  data  as¬ 
sociations  on  the  second  and  third  scans,  subject  to  our  suboptimality  assump¬ 
tions.  Clearly,  only  the  reasonable  hypotheses  result  with  a  small  gating 
parameter  value  which  will  pass  the  gate  of  10.  For  instance,  the  hypothesis 
which  associates  return  1-2-3  (on  the  third  scan)  results  in  1.21.  In  con¬ 
trast,  the  combination  of  1-4-5  is  given  the  value  of  2730.34  which  will  be 
rejected  by  the  gating  logic.  In  this  case,  the  large  parameter  value  is 
effected  by  the  tracking  filter  which  limits  the  "allowable"  distance  between 
consecutive  returns  which  may  imply,  in  turn,  an  unreasonably  high  target 
velocity  or  position  measurement  error. 
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TABLE  6-1. 


AN  EXAMPLE  OF  GATING  EFFECTIVENESS 
(SEE  FIG.  6-18) 


Scan 


Return 


1  1 

2  2,  A  (false  return) 

3  3,  5  (false  return) 


Scan  2 

Associations 

Scan 

1  Scan 

2  Gating  Parameter 

1 

2 

1.05 

1 

A 

A. 60 

Scan  3 

Associations 

Scan  1 

Scan  2 

Scan  3  Gating  Parameter 

1  2  3 

1  2  5 

1  A  3 

1  A  5 

1  Miss  3 

1  Miss  5 

2  3 

2  5 

A  3 

A  5 


1.21 
761A.61 
1236.19 
2730. 3A 
1.01 
1A  .85 
1.19 
AA.AO 
1.59 
39.93 
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SECTION  7 
CONCLUSIONS 

In  this  effort,  an  advanced  algorithm  for  surveillance  radars  has  been 
developed  and  demonstrated.  The  algorithm  can  track  multiple  tracks,  identi¬ 
fying  false  alarms,  accepting  missed  detections  of  targets  in  track,  initia¬ 
ting  tracking  of  new  detected  targets,  and  terminating  tracking  of  targets 
which  are  no  longer  detected.  It  is  based  on  the  statistical  estimation 
theoretic  concept  called  multiple  model  adaptive  estimation  and,  in  the  ab¬ 
sence  of  limitations  on  data  storage  and  processing  rate,  it  is  optimal.  Such 
limitations  forced  an  implementation  which  is  suboptimal  in  theory,  but  which 
closely  approximates  the  performance  of  the  optimal  algorithm  due  to  the 
resiliency  of  the  multiple  model  approach  to  adaptive  estimation. 

In  the  course  of  developing  the  algorithm,  two  models  were  developed  of 
the  multitarget  radar  environment,  one  suitable  for  algorithm  design  and  the 
other  suitable  for  simulating  inputs  to  a  tracking  algorithm.  Both  models 
stand  on  their  own  as  well  as  serving  in  the  development  and  evaluation  of  the 
tracking  algorithm.  The  development  of  the  model  is  firmly  grounded  in  the 
existing  body  of  theoretical  results  on  multiple  model  adaptive  estimation  and 
serves  as  a  good  example  of  the  application  of  these  results  in  its  progres¬ 
sion  through  cases  of  increasingly  complex  environments.  The  reorganization 
of  the  optimal  tracking  algorithm  for  the  most  complex  case  is  particularly 
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significant  since  it  sets  the  stage  for  practically  implementing  the  algorithm 
without,  at  that  state,  compromising  its  theoretical  optimality. 

The  impossibility  of  implementing  the  optimal  algorithm  with  finite  limi¬ 
tation  on  data  storage  and  processing  was  clearly  demonstrated;  for  quite 
reasonable  scenarios,  the  storage  and  processing  requirements  grow  exponen¬ 
tially  or  faster  with  the  number  of  scans.  But  by  starting  from  a  complete 
optimal  algorithm  and  proceeding  systematically  to  develop  rules  for  ignoring, 
pruning,  and  merging  hypotheses,  we  have  begun  to  develop  a  suboptimal  track¬ 
ing  algorithm  which  avoids  such  problems  while  maintaining  performance  which 
is  approximately  optimal.  The  rules  resemble  the  ad  hoc  rules  incorporated 
for  common  sense  reasons  into  many  practical  tracking  algorithms  in  use  today 
and  suggest  how  one  might  profitably  proceed  to  improve  upon  the  presently 
used  rules;  common  sense  does  not. 

The  tracking  algorithm  was  implemented  in  a  flexible  and  general  manner 
to  allow  experimentation  with  various  rules  for  keeping  the  algorithm's  data 
storage  and  processing  requirements  within  reasonable  limits  or  even  to  use 
the  full  implementation  of  the  optimal  algorithm  for  a  limited  number  of 
scans.  A  consequence  was  the  development  of  a  pair  of  data  structures  of 
great  utility  for  organizing  the  data  storage  and  processing  of  a  multitarget 
tracking  algorithm:  the  association  list  and  the  global  hypothesis  list. 

While  not  exploited  in  the  actual  implementation  due  to  the  inherently  serial 
nature  of  the  computer  use^.  (VAX-11/780),  there  are  ample  opportunities  for 
parallel  computation  provided  by  these  data  structures  and  the  procedures 
manipulating  them. 

The  power  of  the  tracking  algorithm  was  well  illustrated  by  the  examples 
of  tracking  performance  provided  in  Section  6,  Target  tracks  can  be  picked 
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out  of  an  environment  so  high  in  false  alarms  that  the  eye  cannot  discern  the 
number  of  targets  appearing  in  10  scans  of  tracker  input  data,  let  alone  their 
tracks.  The  resiliency  of  the  algorithm  is  demonstrated  by  the  Inclusion  of 
the  truth  hypothesis  at  or  near  the  head  of  the  global  hypothesis  list  in 
every  case  despite  their  stressing  natures.  Usually,  the  other  hypotheses  are 
trivially  or  transiently  different  from  the  truth  hypothesis.  Sometimes  a 
false  track  was  hypothesized,  but  these  tracks  are  distinguished  by  being 
based  on  few  and  recent  detections;  potentially  spurious  tracks  can  be  recog¬ 
nized  by  using  that  criterion.  The  only  notable  exception  occurs  with  the 
transposition  of  tracks  in  Fig.  6-17.  In  such  a  congested  environment,  such 
crossovers  are  possible  (e.g.,  a  fighter-on-fighter  engagement)  and  measured 
radar  positions  do  not  by  themselves  provide  the  information  needed  to  recog¬ 
nize  them.  Knowledge  about  typical  multiaircraft  flight  patterns  under  vary¬ 
ing  circumstances  (engagements  versus  patrolling)  are  needed;  such  knowledge 
is  not  readily  captured  in  an  automated  target  tracker. 
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